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ABSTRACT 


The thermal elastohydrodynamic problem of contacts sub- 
jected to rolling and sliding is investigated. The governing con- 
tinuity and momentum equations have been combined to yield a modified 
form of the Reynolds equation. The viscosity and density of the 
lubricant are assumed functions of temperature and pressure. The 
coupled equations of elasticity, Reynolds and energy equations for 
the fluid and the energy equations of the solids are solved iteratively 

U 


for values of the non-dimensional speed parameter Bs ranging from 


ER 
10712 to 10719 and for slip ratios ranging between 0 to 0.25. The 


numerical scheme solves the ‘isothermal’ problem first and the results 
are used as an input to the thermal problem. Comparison with other 
theoretical solutions and experimental results show that the results 

of film thickness obtained in the present work correlate more 
favourably with experimental work. The temperature distribution in 

the lubricant, the surface temperature of solids and the tractive 
forces are also given. The agreement between the calculated and the 
experimental values on drag force available in literature is moderately 


good. 
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CHAPTER I 
INTRODUCTION AND REVIEW OF LITERATURE 


1.1 Statement of the Elastohydrodynamic Problem 

Elastohydrodynamic lubrication may be defined as the study of 
those situations where the elastic deformation of the bounding solids 
has a very significant role in the hydrodynamic lubrication. In most 
machines forces are transmitted from one component to another by means 
of large effective bearing areas. In several applications, as in gears, 
cams and roller bearings, the contact is limited to a point or line. 
The present study pertains to such situations and the geometry envisaged 
is that of two cylinders in contact along a generator. 

In order to fully understand the problem it is necessary to 
consider the features of such contacts. The size of the Hertzian ioe 
zone of contact between two elastic solids is representative of the 
region of effective pressure generation in elastohydrodynamic contacts. 
This is of order 107% inches. The usual surface speeds are of order 
10* to 103 inches per second. The maximum contact pressures are of 
order 10! to 102 tons per square inch. Film thicknesses are of order 
10-5 inches. 

The distinguishing feature of elastohydrodynamic lubrication 
is the elastic deformation of the solids. This is significant compared 
to the lubricant film thickness. The stress distribution determines 
the local shape of the elastic solids. This in turn determines the 


shape of the lubricant film. Thus, the isothermal elastohydrodynamic 


70 vhuds Sa at 
tytioa anthat 
ai i wai 
Sat r vd oh. 


Wy YIAAewoI6 
(1) oulsd7 
ais 1 vise 
»a'oeg005 j3iman 
rab% 10. Sz 
o 8 i 
1 ob 2 he ee 


engines 


vhovhydetents nl coltagegtoy sreeestq syraietts le aol; 


agzassuls tis olis Wet oaupe tsq wid I9OL +08 ‘or 


voeeiiaans Seennvborhuiedia ls lo simwsB44 snide bugalinihh ont 


iadul olmeovhorbyil. ofa al. sks foe DTT hiate vray 8 eet 


ali baw anotiausite tue of 4aladisq VouIe Gpeamag 


9) mutderg ata boedershav ¥biot oF 1shao nt 


mg V3s 709 mir i cee alt  bNO5S2 TOY esionl Ey 67 Or. 




















r 
1 
| 770 ATit) 
mrran eh 70 W1IVES <A POTTING Th : 


mofdoxsY ote rvporbyHiasesla ait Io tosee3are f.£ 


d van néitestodul obeenybosoyderaeds 


a 
tents BAY Brsde asotjeusie seona 


iroish of3 


@ u 
aumes st6 peti besiiaegas? ogs 209307 SaRRe 
— 4 ian 


sotluqe [stevee 67 .eno2b yoties sviesstie sgrel Te 


“+ Hevea! eb ssesaco ofy , ABGESES tIot bas ae 


a 


ia tINSy & or ILG 25837002 iF erabailys ows to. i at 
' 


aos A» tove Yo saves? aa Tet 


43] 
' 


ingan'h Al abtfarn a) s90 le “ows ysawged teaknigod to ee 


~ 


o2 esonltue lRost silt sotinnt S-0¢ Tobt0 to ak ‘ag z 


~ setadioant * 


- 


U 






. 


pina . ace hen oe asa cisunle af2 ® 





problem calls for the simultaneous solution of the continuity and momen- 
tum equations for the fluid and elasticity equations for the solid. 

The very high contact pressures have a marked influence on the 
properties of the lubricant. In particular, viscosity of the lubricant 
increases enormously at high pressures and the compressibility of the 
lubricant is also significant and cannot be ignored. 

When the effect of temperature change is considered, the 
problem becomes extremely difficult and very little theoretical work 
has been done. If the opposing surfaces have nearly equal velocities, 
the problem can be idealised to be an isothermal problem. This has been 
the usual approach. If the velocities of the mating surfaces are widely 
different, viscous dissipation becomes important. Further, temperature 
changes have a significant influence on viscosity. Thus, the basic 
differential equations for a complete study of elastohydrodynamic lubri- 
cation are the continuity equation, the momentum equations, the energy 
equations for the fluid, the energy and elasticity equations for the 
solid and the equations of state for the lubricant. These equations 
are interdependent through the physical properties of the lubricant and 
solid. They are extremely complicated and several simplifying assump- 
tions have to be made before a solution can be obtained. 

Experimental work has revealed the effect of temperature upon 
film thickness to be little. One of the most important conclusions of 
elastohydrodynamic analysis and experiment up to date is that the 
viscosity of the lubricant in the vicinity of the inlet to the contact 
controls the film thickness within the contact. The inlet viscosity of 


the lubricant is determined mainly by the temperature of the surrounding 
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solids. This temperature is not sensitive to the heat generated in the 
fluid within the contact. Thus, film thickness can be calculated 
accurately by neglecting energy equations [2]. Then the effects of 
temperature can be analysed. 

In the present study two rotating elastic cylinders with a 
thin lubricant film inbetween are subjected to a heavy load. Con- 
sidering viscosity and density of the fluid as functions of temperature 
and pressure, and assuming the elastic properties of the solids in con- 
tact, the film thickness, the pressure and temperature distribution 
within the lubricant film and the temperature distribution in the solids 
are computed for different rolling and sliding velocities. The fric- 


tional forces are also calculated under these conditions. 


1.2 Review of Relevant Literature 

The initial interest in elastohydrodynamic lubrication was 
generated by the study of lubrication in gears. Martin [3] made a 
theoretical approach to the problem of lubrication of rigid cylinders. 


For an isoviscous incompressible lubricant he derived an expression of 


the form 
h U; + Uo Uu 
kale Oe aah ins 
R 2245 HU, W 4.9 W e ° ° ° e Ge) 


Substituting representative values for viscosity, load and speed, it is 
found that the thicknesses are too small as compared to actual thick- 
nesses. 

Peppler [4] considered the elastic distortion problem 
assuming an isoviscous lubricant. He came to the conclusion that maxi- 


mum oil pressure cannot exceed the maximum Hertzian pressures. This 
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conclusion is not true for high speeds. Meldahl [5] examined the effect 
of high pressure on film shape and pressure profile for an isoviscous 
lubricant. He derived expressions for the elastic displacement of a 
semi-infinite elastic solid subjected to an arbitrary surface loading. 
He tried to solve the elastic and the hydrodynamic equations together, 
using iterative methods. Convergence was poor and considerable compu- 
tational effort was required for a single solution. Nevertheless, it 
was a move in the right direction. 

Gatcombe [6] took into account the effect of pressure on vis- 
cosity. He used an exponential relation and solved the Reynolds 
equation. His calculated film thicknesses were higher than previous 
calculations. Hersey and Lowdenslager [7] employed a parabolic vis- 
cosity relationship and arrived at results similar to that of CGatcombe 
[6]. Cameron [8] and McEwen [9] employed a boundary condition assuming 
cavitational effects and improved upon the previous workers' approach. 
Block [10] gave a mathematical reasoning proving the existence of a 
minimum film thickness for lubricants obeying an exponential pressure 
viscosity law. Thus far it has been shown that elastic distortion and 
viscosity pressure effects could separately account for modest improve- 
ments in minimum film thickness predictions. 

Grubin and Vinogradova [11] examined the combined effect of 
both elastic distortion and viscosity pressure effects. Grubin made a 
simplifying assumption that the solids adopt the form of dry contact. 
Pressure at entry to the high pressure zone was assumed very high. 
Under these assumptions Grubin calculated the separation of the solids 


within the Hertzian contact zone. This approach eliminated the need 
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of solving the elastic equations since the Hertzian form of dry contact 
was accepted. This assumption is particularly valid at high loads, 
where the hydrodynamic film thickness in the high pressure region is a 
small proportion of the local elastic displacement. The important con- 
tribution of this method was that an approximate formula for film 
thickness for highly loaded contacts could be derived whichwas 2 orders 
of magnitude greater than Martin's [3] expression. Grubin's [11] work 
further postulated the existence of a pressure spike at the outlet end 
of the Hertzian region. 

The search for accurate formulae for pressure distribution 
and film shape has proved very interesting. Weber and Saalfeld [12] 
obtained a closed form solution to the elastohydrodynamic problem for 
constant and pressure dependent viscosity fluids. Their solutions are 
restricted to very small deformations. Their results indicate the 
change of film shape with load. Dowson and Higginson [13] developed an 
interesting numerical method to solve the isothermal elastohydrodynamic 
problem. They found the film shape corresponding to an initially 
assumed pressure distribution by solving the Reynolds equation. For 
the same pressure distribution the elastic film shape was determined by 
solving the elastic equation. The elastic film shape and the hydro- 
dynamic film shapes were compared. If there was sufficiently close 
agreement, the pressure and film profiles were accepted. If not, the 
pressure curve was modified and the procedure repeated till there was 
acceptable agreement. The deficiency in the method was that the adjust- 
ment of the pressure curve had to be done manually and required 


experience and judgement. Further, the results did not reveal any 
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pressure spikes, which was due to the very low velocities for which the 
investigations were carried out. Ina later paper Dowson and Higginson 
[14] investigated the effect of the material properties on the iso- 
thermal elastohydrodynamic situation. They found that an outlet 
pressure peak existed for realistic values of speeds and loads. The 
magnitude of the pressure peak varied very slightly with load but 
markedly with speed and material properties. The centre-line film 
thickness varied hardly with the load and was significantly related to 
the product of speed and inlet viscosity. Archard, Gair and Hirst [15] 
developed an iterative procedure for the isothermal elastohydrodynamic 
problem. The lubricant was assumed to be incompressible. Its viscosity 
was assumed to be dependent on pressure only. They divided the pressure 
region into four regions and employed the inverse hydrodynamic procedure. 
They verified the earlier predictions in a number of cases. The film 
thicknesses they calculated are uniformly higher than those of other 
workers and available experimental values. The method they developed 
is valuable in that the whole calculations can be carried out in a com- 
puter without the need for human intervention which was necessary in 
the methods adopted by Dowson et al. Archard and Kirk [16] examined 
the lubrication of point contacts and came to the conclusion that at 
heavy loads elastohydrodynamic theory had to be applied to achieve 
acceptable results. Dowson, Higginson and Whitaker [17] investigated 
the effects of speed on film thickness. They also took into account 
lubricant compressibility. Speed of rollers had sipnificant 

effect on film thickness and the effect of compressibility of the 


lubricant on film thickness was little. Compressibility had an 
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effect on the pressure peak. Stephenson and Osterle [18] also came to 

similar conclusions. The numerical scheme adopted by these authors is 

however useful only for a narrow range of loads. Dowson and Whitaker 

[19] have also indicated easy methods to determine whether a particular 

problem is to be treated as a rigid cylinder or as an elastohydrodynamic 

problem. Herrebrugh [20] examined the isothermal problem from a mathe- 
matical viewpoint. For an isoviscous lubricant he combined the Reynolds 
and elasticity equations to yield a Fredholm equation of the second 
kind. For a number of loading conditions he has presented numerical 
solutions. In the range covered, he did not find pressure spikes in 

the pressure distribution, showing that spikes occur if the viscosity 

increases rapidly with pressure. 

When temperature variation and its effect on viscosity and 
density of the lubricant is taken into account, the problem becomes 
much more difficult. This is indicated by the limited literature avail- 
able on the subject. Cheng and Sternlicht [21] examined the problem of 
thermal elastohydrodynamic lubrication of rolling and sliding cylinders. 
They made an assumption of an average effective visccsity across the 
film thickness. Instead of hunting for the spike position, they fixed 
the location of the pressure spike and their numerical procedure was to 
find the corresponding rolling velocity. Their important conclusions 
were that: 

(a) Pressure peak existed near the downstream end of the film for all 
heavily loaded cases and position of the spike was dependent upon 
rolling velocity; 

(b) Temperature did not reduce the pressure peak but increased it 


slightly; 
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(c) Temperature had a moderate influence on film thickness. 

The compressible work term used by these authors was in error as was 

pointed out by Dowson [22] in his discussion of the paper. In a sub- 

sequent paper Cheng [23] removed the restriction of average viscosity 
across the film. His form of Reynolds equation has been found to be in 
error by this author. His important conclusions were: 

(a) The position of the peak is strongly influenced by speed, load 
and inlet viscosity; 

(b) No significant difference existed between isothermal film profiles 
and thermal film profiles for various ratios of slip; 

(c) The frictional force was strongly influenced by the temperature 
rise. 

Dowson and Whitaker [24] have presented results for various sliding 

velocities, keeping the rolling velocity constant. In their numerical 

solution, the property expressions of the lubricant are not differen- 
tiated or integrated. This necessitates finding thermal coefficients 
in order to use the modified Reynolds equation. Their main conclusions 
are: 

(a) The effect of temperature on the film thickness is negligible; 

(b) Sliding speed has little effect on the film thickness; 

(c) Height of the pressure spike is reduced by sliding and its location 
moves towards the centre of the Hertzian zone. The spike becomes 
more gentle and rounded. 

One important common feature of the solutions offered by the 

foregoing three references [21, 23, 24] is that all these assumed a 


surface temperature equation of the solids as a boundary condition on 
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temperature for the fluid. The temperature distributions within the 
solids were not investigated. When it is realised that experimental 
measurements of temperature have necessarily to be done beneath the 
surface, the temperature distribution in the solids assumes importance. 
Shear stress is significantly dependent on temperature. 

O'Donoghue and Cameron [25] correlated all available data on 
friction between hardened steel discs and it is possible to see the 
effect of solid temperature on friction. In another paper these 
authors [26] investigated the mechanism of scuffing between gears. They 
came to the conclusion that scuffing takes place when the surface 
attained the transition temperature of the particular metal-oil-metal 
combination. If this temperature is reached then scuffing will ensue 
(even if the surfaces are fully lubricated) due to chance contacts. 

Blok [27], Carslaw and Jaeger [28], Jaeger [29], have studied 
the theoretical problem of a moving heat source of constant strength. 
It is usual in elastohydrodynamic studies to consider the source 
strength over the Hertzian width as constant. The question of partition 
of heat between two bodies has been studied by Allen [30] and Cameron, 
Gordon and Symm [31]. These authors took as their system a source of 
uniform strength and put in the condition that at all points in the 
zone of contact the surface temperatures of the solids must be equal, 
They have derived expressions for temperatures in terms of the source 
strength. Francis [32] has derived an analytic expression for the 
interfacial temperature in a sliding circular Hertzian contact, where 
one surface was stationary and the other moving. He has introduced the 


method of the harmonic mean for matching the interfacial temperature. 
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It has been suggested that the non-newtonian behaviour of the 
lubricant might influence the film thickness in elastohydrodynamic con- 
tacts. Cheng and Orcutt [33] considered a visco-elastic model of 
Maxwell and found that the shape of the pressure profile obtained was 
closer to the experimental results. Milne [34], Tanner [35] and Bell 
[36] have come to the conclusion that the load capacity of a Maxwell 
fluid was less than that of a Newtonian fluid. Bell [37] derived a 
formula for film thickness for a Ree-Eyring fluid. Chow and Saibel [38] 
analysed the isothermal problem of a heavily loaded line contact of 
rollers in the presence of a Maxwell lubricant. They have found that 
the non-Newtonian effect had the tendency to flatten the contact region. 
The load capacity was reduced. The shape of the pressure curve did not 
vary significantly. The existence of a spike was found as in the case 
of newtonian fluids. 

A number of experimental investigations have been made for 
the problem under study. A brief review of these is presented to indi- 
cate the extent of the correlation between available experimental re- 
sults and theoretical findings. Merrit [39] was the first to design 
and build a disc machine which simulated gear tooth contact conditions. 
His measurement of the coefficient of friction under various sliding 
conditions presented valuable design data. Lane and Hughes [40] studied 
the oil film formation between gear teeth using electrical resistance 
method. They found that under sliding conditions film thickness was 
less than under pure rolling conditions. Crook [41] investigated the 
possibility of a quantitative evaluation of the film thickness by a 


knowledge of the specific resistance of the oil. He found that even 
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very small quantities of moisture vitiated the results considerably. 

El Sisi and Shawki [42] tried to overcome the effect of moisture by 
adding 4 per cent sodium petroleum sulphonate. But since they used 
rather heavy currents, their results were one order of magnitude higher. 
Leach and Kelley [43] have used electrical resistance methods across 
the contact zone to identify lubricant failure point and the effect of 
deposit forming additives. They have come to the conclusion that the 
load capacity of a lubricant varies inversely with reference to speci- 
men temperature and that lubricant failure for any lubricant-material 
combination occurred at a constant critical temperature. This tempera- 
ture did not depend on the range of load, sliding or rolling velocities 
or film thickness. Tallian, Chik, Huttenlocher, Kamenshine, Sibley and 
Sindlinger [44] have found that there is a good order of magnitude 
agreement between measured and theoretical values of film thickness. 
Significant wear occurred only in regions where the film was interrupted. 
All the experimental investigations using resistance method have quali- 
tatively confirmed predictions of the elastohydrodynamic theory. 

Sibley and Orcutt [45] devised a method which consisted of 
directing a mono-chromatic beam of X-rays tangentially at the contact 
of the lubricated discs. The amount of the X-ray beam passing through 
the gap was measured by a Geiger counter and was used to determine film 
thickness. The wave length of the X-rays was so chosen that the lubri- 
cants were quite readily penetrated but the steel surfaces absorbed the 
beam entirely. The elastohydrodynamic theory fitted the experimantal 
results well up to an oil film thickness of 10~° inches. At lower 


thicknesses (the lowest thickness measured was 3 micro-inches) it 
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deviated twenty per cent below theoretical. The authors attributed this 
deviation to the non-Newtonian behaviour of the lubricant at high shear 
rates. The X-ray method relies on the direct transmission of X-rays 
through the contact which limits it to discs. 

Measurement of the electrical capacity between two surfaces 
separated by an oil film has been found to give good results. Crook 
[46] used a modified disc machine so that the oil coming from the con- 
tact zone remained on the surfaces and was carried under a flat plate. 
The capacitance between each plate and its disc was measured. The use 
of the subsidiary plate eliminated uncertainties of the geometry of the 
contact zone. From the capacitances the rate of oil flow was deduced 
and thus the thickness of the film could be calculated. The method is 
applicable for both rolling as well as rolling and sliding conditions. 
The important findings of these measurements were that film thickness 
at high loads varied little with load, less than proportionately with 
speed and greatly with the temperature of the surfaces. The tempera- 
ture of the oil in the high pressure region had negligible influence on 
the film thickness. Crook thus came to the conclusion that film thick- 
ness is largely determined by the conditions on the entry side of the 
conjunction ahead of the region in which the viscous losses and heating 
of the oil became appreciable. In a second paper Crook [47] investi- 
gated the effect of viscosity and speed on film thickness. He found 
that the viscosity of greatest importance with respect to film thick- 
ness was the viscosity of the oil at the surface temperature of the 
discs. It was also found to depend upon the mean peripheral speed of 


the discs. This was found to be valid even when sliding conditions 
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prevailed. Ina third paper Crook [48] found that under sliding con- 
ditions conduction to the solid surfaces to be the important mechanism 
of dissipation of heat. He has also developed expressions for frictio- 
nal traction. In a fourth paper Crook [49]: has described the measure- 
ment of friction in disc machines. He found that rolling traction was 
independent of load and was proportional to the thickness of the film, 
Sliding traction was dependent upon sliding speed. Friction was found 
to- increase with sliding speed to a maximum and then full. The four 
papers [46, 47, 48, 49] have been described in some detail because they 
represent a landmark in the experimental investigation of elastohydro- 
dynamic lubrication in line contacts. Dyson and Wilson [50] have 
investigated the effect of high slide/roll ratios on film thickness, 
using electrical capacitance method. They found that as the slide/roll 
ratio increased above a certain level, experimental values of the film 
thickness became greater than the values predicted by isothermal theory. 
They explain this on the principle that there is a temperature variation 
of oil across the film at the inlet zone and this becomes important at 
high sliding rates. 

Kannel, Bell and Allen [51] have reported two different 
methods for measuring pressure distributions in rolling contact. The 
first method was to use an X-ray technique and pressures were inferred 
from observed deformations of the discs. The second method used a strip 
of manganin as a pressure transducer to measure pressures in the con- 
tact. The important feature of the measured profiles were that they 
did not show the pressure spike predicted by theory. The authors thenm- 


selves have posed the question whether lack of accuracy could have been 
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a cause of the absence of pressure spike. Qualitatively the pressure 
profiles have agreed with the theoretical predictions. Kannel [52] 
refined the manganin pressure transducer technique so that pressures 
could be measured between steel discs. The measured pressures showed 
general trends similar to those predicted by theory. For heavily 
loaded rolling-contact conditions, a slight pressure spike has been 
reported. Longfeld [53] used oil pressure tappings to measure directly 
the pressure in the contact. Facilities were available enabling 
measurements of pressure over the entire contact area. At greater 
loads, rudimentary peaks were observed on the downstream end of the 
contact. Orcutt [54] has reported a detailed experimental study of the 
conditions occurring in the conjunctive region of two lubricated 
cylindrical discs which roll or roll and slide on their peripheral 
surfaces. Platinum transducers were used to measure temperature. 
Capacitance method was adopted to measure film thickness. Manganin 
pressure transducers measured the pressure. Experimental pressure re- 
sults indicated that significant pressures were not generated until 
about two Hertzian contact zone half widths ahead of the line of centres. 
Experimental pressure profiles were more rounded without the sharp 
pressure peak which may be due to the rheological effects. The 
measured value of the film thicknesses was about 60 to 70 per cent of 
theoretical values. The shapes of the deformation profiles differ from 
the calculated profiles which are flatter than the measured ones. The 
measured temperature rise is less than the calculated values. However, 
these measurements indicate qualitative agreement with theoretical re- 


sults. 
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Optical interference technique has been used for the measure- 
ment of film thicknesses. The advantage of this method is that the 
measurement is made at the actual contact and that it is independent of 
external calibration. The difficulties are two-fold. Firstly, most 
transparent materials such as glass have refractive indices very close 
to that of normal lubricants. Secondly, glass is a bad bearing material. 
If perspex or some other plastic is used then the modulus of elasticity 
and yield stress are small so that the pressures generated are too low 
for the pressure-viscosity effect to occur. Cameron and Cohar (554 
employed glass with a high refractive index which enabled the colours 
of the fringes to be photographed. The opposing surface was a one inch 
steel ball. The pressures were high enough to get an elastohydrodynamic 
lubrication. The experimental results gave the exit constriction pre- 
dicted by theory. Very recently, Westlake and Cameron [56] have used 
optical elastohydrodynamic techniques for testing a wide range of lubri- 
cants. The pressure viscosity coefficients agree well with those de- 
termined in conventional high pressure viscometers. The effect of aging 
on the oil could be studied. Sanborn and Winer [57] have used the 
optical interference technique for the study of the rheological effects 
in sliding elastohydrodynamic contacts. They used a steel sphere and 
synthetic sapphire in their apparatus. Film thickness profile was not 
affected by rapid application of the load. These authors also measured 
the temperature at the inlet and it was found that there was only a very 
slight increase in temperature after a long time which affected the film 
thickness only slightly. In a companion paper Sanborn and Winer [58] 


have described traction measurements. The traction values obtained 
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were primarily a function of the Sliding velocity. Large variations 
in fluid composition and inlet viscosity had little influence on the 
traction force. Rapid application of the normal load also had negli- 


gible effect on the traction force. 


1.3. Approach to the Thermal Elastohydrodynamic Problem 

Combining the continuity and the momentum equations for a 
compressible fluid, a modified form of Reynolds equation has been 
derived in Chapter II. The energy equations for the fluid and solids 
have also been detailed. The elastic equation and the applicable 
Hertzian formulae have been listed in Chapter II. The problem under 
study was solved in two steps. The first step was to idealise the 
problem so that the fluid properties were dependent on pressure only. 
This is usually known in the technical literature as the 'isothermal' 
problem. The manner in which the coupled elastic and Reynolds equations 
have been solved has been shown in detail is Capcer =! Di = ihisestep 
yielded the pressure distribution and film thickness. 

Using the well documented assumption that the film thickness 
was not sensitive to the temperature of the lubricant in the high 
Pressure region, the second step was to devise an iterative numerical 
procedure which took into account the variation of the fluid properties 
with reference to temperature as well. An implicit finite difference 
technique has been developed to solve the energy equations of solid- 
fluid-solid together. These numerical procedures are the subject matter 
of Chapter IV. 


The solution to the momentum and energy equations enable the 
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calculation of shear stress and the tractive forces. The pertinent 
expressions have been listed in Chapter IV. Results for different 
values of rolling velocity for a constant load and variable slip ratios 
have been displayed. The rolling velocities varied from 30 inches per 
second (slow) to 400 inches per second (moderately high). Discussion 
of these results and conclusions thereof form the subject content of 
Chapter V. 

While some of the assumptions have been indicated at the 
appropriate places, it has been considered pertinent to list the general 
assumptions usually made in the derivation of governing equations. 


This is done in the next article. 


1.4. General Assumptions 

For the contacts under consideration in the present study, 
the loads are transmitted through lubricant films of very short length. 
For this reason the undeformed solids can be adequately represented by 
cylinders in the region of the contact zone. Following Dowson and 
Higginson [59] it is possible to replace the two cylinders by a geo- 
metrically equivalent cylinder near a plane. 

The following assumptions are made in the development of the 
governing equations. 
1. Body forces are neglected. 
2. Pressure is constant across the film. 
3. Inertia forces are neglected. 
4. Curvature effects are small and can be neglected. 


5. The length of the contact region is orders of magnitude smaller 
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than the radii of the rollers. Thus, surface deformation of the 
rollers can be approximated by those of semi-infinite solids sub- 
jected to the same normal load. Elastic deformations due to sur- 
face shear is neglected. 

Viscosity is a function of pressure and temperature only and not 
dependent on shear rate. 

Side leakage is neglected. 

There is no slip at the boundaries. 

Flow is considered laminar. 

The thickness of the lubricant film is orders of magnitude smaller 
than the length of the contact. 


The problem is two-dimensional and steady-state has been reached. 


7 
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CHAPTER II 


FLOW GOVERNING EQUATIONS 


2.1. Modified Reynolds Equation 

One of the major assumptions in the two-dimensional Reynolds 
equation, commonly used, is that viscosity is constant throughout the 
thickness of the film. Cheng and Sternlicht [21] improved their solu- 
tion by using an amended form of the Reynolds equation where an 
"effective' viscosity of the lubricant considered constant across the 
film was used. Cheng [23] derived a modified Reynolds equation which 
removed the assumption of the 'effective' viscosity. Unfortunately, 
certain errors have occurred in Cheng's analysis and these seem to 
have been left unnoticed so far. In this article, the thermal form 
of the Reynolds equation is derived and compared with Cheng's [23] 
expression. To facilitate direct comparison, Cheng's notation has been 
adopted. 


Neglecting inertia, the momentum equation is 


ope 22, oe 

dx By “u oy (251) 
The boundary conditions are 

u(x,o) = Uj) u(x,h) = U9 (22) 

p(-») = p(x*) = 0 = Bex) (2.3) 
Integrating equation (2.1) once with respect to y 

ou. op ¥ j AG) (2.4) 


dy dxu u 
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Integrating equation (2.4) once again with respect to y 


u= 


a) 


y y 
{ e dy + A(x) f . dy + B(x) (2.5) 
0 O 


Using boundary condition u(x,o) 


U,; in equation (2.5) one 


obtains 


B= Uj (2.6) 


Using boundary condition u(x,h) Uy in equation (2.5) one 


obtains 


U2 = gel Lay+am f tay+u, 257) 


1 ly i 
ula os sok = dy (2. 8a) 
ue aA 
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1 2 
aie ae. ae (2. 8b) 
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a = mae I : dy (2. 8c) 
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From (2.7) 
h2 
a Na lair 10 
F e 
A(x) = “he se (2.9) 
We 
We define an effective density p = 0 (x) C2) 
such that 
ae + 
mass flow = Q = 60 f u dy (2501) 
fe) 
Integrating by parts equation (2511) 
» h h Ju 
Oe = Ton anuy. oa Yio aye 
= dp /” y? Ny 
=p {Ugh - <2 f mpd yann A es) 
O fo) 
hy 
ae feed Dene te dp bh? 
me He {U5h 3u" dx a 2 (Uj-U}- dx at} 
e e e 
= d hou, h3 yee 
=p {Uoh ee Orns - 3? = ome ih bie 
e e e 
f tidp 13 toh 
=o {Ush - Tide 4, = SU5h s+ = (Uj+U}) } 
US het) 
Ge pele Gp ie cee 


Dowson and Wiitaker [24] ignored variation of density across the film 


thickness. 
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Following Cheng and Sternlicht [21] and Dowson and Whitaker [59] 


we adopt the boundary condition that 


Is 


at x = x*_ p 














Then 
eet 

Q ae oy Nae eead Sll sh (2a3) 
Equating mass flow 

- U,+U> U,t+U> het} 

= = Lid 

* * *) ee +. GP 
eee) Bey a Net GS, 
h*p* 
6u_ (U,+Uo) {to - —— to*} 
PE a seas deme tiD =e 2 Oy 
dx t7h 


The expression (2.14) is to be contrasted with Cheng's [23] 


expression which reads 





6u (U,+U>) sya o> 
ONS ee pe (2.15) 
dx tyh ho 


This is true only if t»)* = 1 and from equation (2.8e) it 


follows that for ty* to be equal to unity 
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Thus, it seems to be tacitly assumed that the temperature at x = x* 
is constant. The validity of this assumption is certainly questionable. 


However, it can be used to start an iterative process in the 


numerical procedure. 


2.2. 'Isothermal' Expressions 
The expressions for the case where lubricant properties are 
only pressure dependent can be easily deduced from equation (2.13). 


For the ‘isothermal' case 





6u (U,+U>) {1 - 
dx h2 


6u (U,4+U>) 


b: Hee ee a4 (2.16) 


2.3. Energy Equations for the Fluid and Solids 
Following Dowson and Whitaker [24] the reduced form of 


energy equation for the fluid is 


oT L SP OU) 2 oe DAG 
ae aie euT ae + Gy +k aya (Zea) 
Boundary conditions are 
T (-~, 0) = T (2.18a) 
oT} 
wat 


ay 'y=0 ~ “1 By aeean (2.18b) 
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The temperature distribution in the solids is governed by [60] 


97T, 977) oT] 
ky Ge + mayo = P16, U1 aot (2.19) 
97T>.  92T> dT 
ko Grea + ty? = Pap U2 Ey (2520) 
Non-dimensionalising the equations as shown below 
eevee ete x oe 3 eon) 
T L h 
s 1 r 
we obtain 
h.? 327, 092T, nh? pee 1 97, 
2 a er ee (2.22) 
inte tbe ay Lic iia tat x 
nL alba se clay 1h MeN Mee 
= = (3) 
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substituting (2.21) in (2.18a) yields 


k h 
pay i al 
r 
usually mie san 


= 0 (10-3) in an elastohydrodynamic contact 
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hence 


hee 
1g ~_ = 
1,2 = © |(200 x 107%)2| = 9 {10 -2| 


consequently in equations (2.19) and (2.20) we may neglect second order 


derivatives with respect to x and obtain 


927) oT} 

ky aye) = Pica! Ee (2.24) 
a-T> aT > 

ko Boys = pares ea (225) 


O'Donoghue and Cameron [26] have obtained very satisfactory 
correlations between experimentally obtained values of temperature and 


calculated values using relations similar to (2.24). 


2.4. Elasticity Equation and Useful Hertzian Relations 


HERTZIAN PRESSURE R 
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For the co-ordinate system shown in Fig. 1, the vertical dis- 
placement of a two-dimensional elastic half space under a distributed 


normal load p(&) can be determined [61] from the following equation 


h=h, +o PCE) In ee =| ai (2.26) 


For the following Hertzian relations [60] 
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1 pie 
a= {( E ‘ie Ep ) = WR} (22) 
Daas = (2.28) 
x2, % x 
Be=apee (l= 32) il) (2.29a) 
= 0 (S[>D (2.29b) 


equation (2.26) yields [15] 


2 +5 
h - h, -i- = (3 (esa fan [al + Ge =) ] 
(JS|>) (2. 30a) 
=) (|= | <1) (2.30b) 


2.5. Property Relations 
Viscosity and density of the fluid are functions of pressure 


and temperature. Following Cheng [23] the following property relations 
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are assumed 





2 2 eee: 
Hae (he expy (epee a ter) CAeub) 
C,P 
Peron {1 + 14, + D.(I-T_)} (ey) 


For the ‘isothermal’ case these relations reduce to 


UW = LS exp (ap) C2135) 
C,P 
anes ue * Tce’ (2.34) 
where a= (a pee (2535) 
s 


2.6. Non-dimensional Procedures 
The equations can be non-dimensionalised by introducing the following 


non-dimensional variables. Within the fluid: 
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Within the solid: 


= y = 
EO =e (2.37) 
a Yr. 


Then equations (2.8a) to (2.8e) are non-dimensionalised. to 
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respectively. 


The equations (2.14), (2.16) and (2.3) are non-dimensionalised to 
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The elastic equation (2.26) is non-dimensionalised to 
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The energy equation for the fluid (2.17) is non-dimensionalised to 
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iT (2.48) 


CQ 
i] 


Joule's constant 


One may also define a modified Reynolds number 


2 
ee fe) 
Re* = (2.49) 
uot 
Equation (2.46) could then be rewritten as 
= By hoe --= dp Wes 1 9¢T 
PrRe* pu — = PrM —— euT — + PrE = (4) += —— (2.50) 
Ox L2 dx he Boy h? oy? 
Boundary conditions given by equations (2.18a), (2.18b) and 
(2.18c) are non-dimensionalised to 
T (-~, 0) = 1 (7 Bia) 
= ke oT 
ie = (—) —|_ (2. 51b) 
dy y = 0 Win ye 9 
= Ke Par 
cliy | (2.51c¢) 
hy Ae ak oy, y=l1 


respectively. 


The energy equations for the solids (2.19) and (2.20) are non-dimensionalised 
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= RS, —— (253) 


respectively, where 


ip ry 
RS, = Siecicy a (254a) 
Lp ,C,, Ube . 
2 2 
RS, = (———_) (— ) (2545) 
K, 1h 


and they denote a reference solid non-dimensional number pertaining to 
bottom and top solids respectively. The boundary conditions on 


temperature are given by equations (2.5la), (2.51b), (2.51c) and by 


igs (Es h,? S uy, (x Bo Si (2355) 


The property expressions (2.31) and (2.32) are non-dimensionalised to 
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respectively. 
For the 'isothermal case' equations (2.33) and (2.34) are 


similarly non-dimensionalised to 
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respectively, where 


(2.60) 


The non-dimensionalisation adopted transfers the original solid-fluid- 


solid field into a rectangular field shown in FiO 2 
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CUAPTE Rael. 


SOLUTION TO THE ISOTHERMAL PROBLEM 


3.1. General Approach 

The 'isothermal' problem requires the simultaneous solution 
of the Reynolds equation (2.40) along with the elasticity equation 
(2.43) taking into account the variation of viscosity and density with 
pressure as given by equations (2.58) and (2.59) respectively. As a 
closed form solution of the coupled integro-differential system of 
equations could not be found, a suitable numerical method had to be 
adopted. 

In relatively low load cases, a straightforward converging 
iterative procedure of assuming a pressure distribution, calculating 
the film shape through the elastic equations and re-evaluating the 
pressure distribution by integrating the Reynolds equation can be de- 
vised [12, 18]. However, for high loads, this procedure fails to con- 
verge. 

The next important consideration is whether to keep a fixed 
centre-~line film thickness and calculate the pressure distribution 
suitable for this assumed thickness or else to assume a given load and 
calculate the centre film thickness. Osterle and Stephenson [18] 
assumed a fixed centre-line film thickness and their results show poor 
convergence even for light loads. Further, even from the design point 


of view, only the load to be transmitted and surface speeds are usually 
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known apriori. In the present work, load is considered known and the 
numerical procedure developed yields the film profile and pressure 
distribution. Earlier work [15, 19] has also revealed that the con- 
vergence is much better when load is taken as an input.< 

It is convenient to divide the contact region into sub- 
regions A, B, C and D, as is shown in Fig. 3. The methods adopted in 
each of these regions to solve the elastic and Reynolds equations to- 


gether have been detailed in the subsequent articles. 


3.2. Initial Considerations for Integration in Region A 


From the consideration that 


SB = Gp ee aD) 


and using the viscosity equation (2.58) we can rewrite equation (2.40) 


in the form 





se gs hy es I i ae (3.2) 


Thus, if in the entrance region A (Fig. 3) one could calcu- 
late h, then for any assumed value of ls equation (3.2) can be 
numerically integrated. For the initial cycle of calculations the film 
profile in region A was assumed to correspond to the Hertzian profile 


of dry contact. This assumption has been adopted by many investigators 
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bine 3S lmandahaci heen used by Grubin and Vinogradova [11] to con- 
Siderable advantage when they derived their well known film thickness 
formula. The applicable Hertzian relations have been listed under 
article 2.3.equations (2.27) to (2.30b). 

A Runge-Kutta fourth order integration technique in double 
precision arithmetic (detailed in appendix A) was adopted to carry out 
the numerical integration. The choice of this method depended on the 
considerations that the method is self-starting and has a low truncation 
error [62] and at any stage the step size can be varied. For numerical 
purposes -@ was located at x = -5a (X = -2.5). This assumption has 
been made on the basis of experimental evidence [525 53i6 Various 
authors have chosen basically similar criteria to fix the value of the 
pressure at the boundary AB. Dowson and Whitaker [19] calculated a 
number of pressure curves, each one corresponding to an assumed ahs and 
selected that pressure curve which can be extrapolated to run smoothly 
into the elliptical pressure profile of dry contact. This called for 
human intervention and judgement. In the present work, following 
Archard, Gair and Hirst [15] it was postulated that had ho been guessed 
correctly in the first instance, the extrapolated value of the pressure 
at x = -0.4 will be the same as the value of the Hertzian pressure at 
that point. At any rate, this is just a first step to get the procedure 
started and the iterative procedure takes care of the shape of the final 
curve. As will be seen later, the values of the final film thickness 
calculated agree quite well with film thicknesses calculated with 
these criteria. Several guesses of Lh may be necessary to meet the 


aforementionec initial criterion but this is easily programmed in a 
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computer since increasing h, decreases the pressures in region A and 
vice-versa. LE be is chosen too small one obtains an error message 
from the computer that logarithm of a negative quantity is attempted 


to be computed... 


3.3. Numerical Integration of the Elastic Equation 

The next step is to calculate the displacements in region B 
due to pressures in region A through equation (2.43). However, equation 
(2.43) has a singularity at & = €. Various authors [P3 de Leis have 
used different devices to remove this etredlarity: In effect, these 
consist of dividing the pressure curve into small segments and consider 
p(—) as a second order polynomial [13] or a first order polynomial 
[18, 15] and formally integrate the assumed pressure function. Wernick 
[63] found that such methods introduced a ripple in the film profile 
which he attributed to the neglect of the effect of curvature in the 
pressure that was not fully taken care of in the foregoing methods. 
His method has been indicated briefly in appendix B. Since there was 
a rapid change of curvature in the pressure profile in the entrance 
region and as will be seen later in the region near the spikes, this 
method was considered the most appropriate and was used. In order to 
have a check on the accuracy of the calculations to be expected of this 
method, the values of the film thickness as calculated through equation 
(2.30a) and through Wernick's method for a Hertzian dry contact are 
compared at discrete points in Table I. 

The displacements in region B due to the pressures in 


regicn A are added te the Hertzian displacements in B. 
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COMPARISON OF ELASTIC DISPLACEMENTS FOR U = 400 and W = 0.32 TON/IN. 
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3.4. Centre-Line Shift Technique 

It will be seen that while the Hertzian pressure profile 
for a dry contact is symmetrical about the centre line, the hydro- 
dynamic pressure profiles are not. For this reason, when the surface 
displacements are calculated taking into account the hydrodynamic 
pressures in region A, the elastic displacements will be tilted with 
respect to the line of centres of cylinders. This leads to an oil- 
film thickness of decreasing height in the direction of surface motion. 
But a feature of the elastohydrodynamic contact is that h must be 
sensibly constant over a major portion of the Hertzian zone and parti- 
cularly in the central region. This problem is conveniently overcome 
by employing what is ome as the centre-line shift technique [64,59]. 


This technique has been indicated in appendix D. 


3.5. Inverse Hydrodynamic Relation 

In region B (Fig. 3), the pressures and viscosities are 
high. A straightforward integration of Reynolds equation in this 
region leads to inaccuracies [19]. It has been found [13, 15, 17] 
better to compare the film shapes calculated by the elastic and hydro- 
dynamic relations and modify the pressure curve to ensure close 
agreement. For this purpose it is convenient to rewrite equation 


(2.40) in the form of a cubic equation 
- . —~ xox 
Foe oe nee eens Gn 
dx 0 


A rapid and accurate solution to the equation could be had by adopting 


a Newton-Raphson iterative technique. This procedure has been 
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described in appendix C. For the initial cycle of calculation 
pressure is assumed Hertzian (except at boundary AB, where the value 
already calculated is used) and the pressure gradients are calculated 
at discrete points in region B. The values so calculated are used in 
the solution of equation (3.3) which yields h at the points chosen. 
The values of h are compared to the elastic displacements at the very 
Same points as calculated through methods developed under article 3.3. 
The differences in the film thicknesses as calculated 
through the elastic and inverse hydrodynamic relations are a measure 
for the amendment of pressures needed at these points. An inverted 
Wernick method was developed to calculate the change of pressures 
needed and this procedure has been indicated in appendix B. With 
these new pressures, pressure gradients can be calculated and the 
cycle starts once again with the solution of equation (5.3)0 Leis 
illustrative at this stage to present the block diagram (Fig. 4) to 
indicate the procedure adopted for the iterative solution of the 
elastic and Reynolds equations in Fegions Avand Bvof Fag. 3.0  Ttewill 
be noticed that one is solving for the centre film thickness by con- 
sidering pressures and displacements in regions A and B and Hertzian 
pressures in C and D. The mathematical argument for the procedure 


adopted is given below. Expression (2.16) can be rewritten as 


hop 
ae (3.4) 





ney dpe s 6U (h - 
u dx 


Differentiating both sides with respect to x and rearranging 
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Input W, U,, Us, FIGURE 4 
E 
ASSUME HERTZ DISTRIBUTION Calculation of hy and integration 


No 
ITER= 1 NCOU=1 


Calculate displacement in A, B, 
C, D using Wernick method. 


procedures for regions A and B of Fig. 3. 
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2 
Fb ae fe ra 
dx Lbegeacl 36 


= 6U Pode ap 
smiths (345) 
‘ ; : : dp i 
For Hertzian pressure distribution a Ofatox =)0° and 
expression (3.5) reduces to 
3 
dh Oden 
antes. 5 | (3.6) 
dx x 0 6Uu, dx“ x = 0 
For Hertzian pressure distribution from expression (2.29b) 
ap ee 
dx2'x = 0 a? 
Thus, 
3) 
dh melee 
dx x = 0 6Uy a? Sl, 
Hence, 1£ h, is small and > dis large (as is the case) 
ey at) 
dx x = 0 


This shows that under the conditions of our present problem, one may 


solve for the regions A and B entirely before going on to the regions 
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3.6. Integration Procedure for Regions C and D 

Different procedures have been Suggested by various authors 
[15, 19, 24] for the integration in regions C and D. The choice of 
the technique depends on the assumption of the existence of the 
pressure 'spike'. Investigating the problem for low loads, Stephenson 
and Osterle [18] did not find any spikes. Herrebrugh [20] investi- 
gating the problem for an isoviscous lubricant also did not find any 
"spikes'. The existence of the spike has however been authenticated 
by other authors working under less restrictive assumptions [15, 19, 
21, 24, 38] and also verified by experiment [52, 53, 54]. Cheng and 
Sternlicht [21] avoided the problem of hunting for the location of 
the 'spike'. They fixed the location and found the corresponding 
rolling velocity which will yield the spike at the chosen location. 
However, in any design problem it is the surface speeds that are known 
apriori. 

In this work, a simpler and quicker method to locate fairly 
closely the approximate location of the spike (within 4% of the semi- 
Hertzian width) was found and this was subsequently refined by an 
iterative procedure. Initially the floating boundary CD was located 
say at x = 0.4. In the regions C and D the pressures are assumed 
Hertzian and zero respectively. For the already found pressure 
distributions in regions A and B and the presently assumed pressure 
distribution in regions C and D the elastic equation (2.43) is used 
to obtain the displacements at discreet points in D. For the dis- 
placements found for discreet points in D, equation (2.40) is inte- 


grated. The integration starts from x = 0.5 with the boundary 
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condition p(0.5) = 0. Simpson's rule is found sufficient for inte- 
gration in this region. Integration is carried out till CDeis: reached. 
The value of the pressure at CD serves as a sensitive indicator for 
the location of the spike. If, for example, CD is to the right of 
the actual location of the Spike, pressure found by the integration 

ae CD is small, If CD is to the left of the actual location, the 
pressure will attempt to pass through an infinite value with an error 
message from the computer that the logarithm of a negative number is 
attempted to be computed. Thus, the direction in which CD has to be 
moved becomes known and CD can be loacted. Pressures in region D 

are calculated for this tentative location of CD. With the present 
pressure distribution in D, the already calculated pressure distri- 
bution in A and B and assumed Hertzian pressure distribution in Ce 

the elastic displacements at discrete points in C are calculated 
through equation (2.43). For the Hertzian pressure distribution in 

C (Fig. 3) solution of the cubic relation (3.3) will yield the hydro- 
dynamic thicknesses at the chosen discrete points in C. As before, 
the differences in the film thicknesses calculated through the elastic 
and inverse hydrodynamic relations are a measure for the amendment of 
pressures needed at these points. The inverted Wernick method 
(described already) is used to calculate the pressure amendments 
needed. It is illustrative at this stage to present the block diagram 
(Fig. 5) to indicate the procedure adopted for the iterative solution 


of the elastic and Reynolds relations in regions C and D. 
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FIGURE 5 


ahaccetocalinn of GD Location of the pressure "spike' 
and the integration procedures 


for regions C and D of Fig. 3. 
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As is shown in the block diagram, the new pressures in C 
are used to calculate the elastic displacements in D. And when these 
displacements are used to integrate for pressures in region D, it may 
become necessary to shift the floating boundary CD to achieve closer 
agreement of pressures. In actual calculations it was found that only 
a very slight shift to the left was necessary to obtain satisfactory 
agreement and between two to four iterations were enough to locate 


the position of the spike within 171000 of the Hertzian semi-width. 


3.7. Discussion on the Numerical Scheme 

One of the contributions of the present study is that a 
completely automatic scheme has been developed which can be extended 
to the study of the thermal system as will be seen in Chapter IV. We 
had to make a few assumptions in developing the scheme and the purpose 
of this article is to discuss the validity of the assumptions made in 
the light of experimental findings and theoretical work. The first 
assumption that was made was that the integration can be started in 
region A from x = -2.5 (x/a = -5.0). The results that are presented 
in the next article show that this is adequate, as the pressure becomes 
significant only from x = -1.0. This is also in agreement with the 
experimental findings of Kannel [52] and Longfeld [53]. Im starting 
the iterative process for finding the centre-line film thickness, it 
was assumed that the extrapolated pressure should match the Hertzian 
pressure at x = -0.4. This consideration is based upon the following 
reasoning. The hydrodynamic equation (3.2) shows that h - ho has to 


be small when the pressures are high. From the elastic standpoint, 
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the Hertzian pressure distribution yields the flatness that is observed 
and verified. Hence the system should be closer to the Hertzian 
Pressure distribution. Since the hydrodynamic equations reveal a 
pressure build-up in the entrance region, a certain amount of pressure 
must be subtracted from the Hertzian distribution to compensate for the 
external excess pressure. The true pressure curve may be expected to 
cut the Hertzian curve at a value of p well below the maximum since 
the depression caused by a force diminishes Slowly with the distance 
and the flatness at the high pressure region has to be maintained to 
correspond with the experimental findings [46, 47, 52, 53]. Thus, for 
the starting of the schemes it was assumed that the extrapolated 
pressure should match the Hertzian pressure at x = -0.4 where the 
pressures are well below the maximum. The iterative procedure improves 
the approximation. 

Different schemes are adopted in regions A, B, C and D to 
find iterative solutions to the elastic and hydrodynamic equations. 
Previous studies by others have shown the necessity for such a step. The 
numerical method to solve for the inverse hvdrodynamic equation, adopted 
in the present work, is new to the elastohydrodynamic studies and has 
been found to be quick and accurate. The inversion of the Wernick method 
to find the amendment to the pressure curve is also a new feature. As 
this method takes into account both the pressure gradients and pressures, 
the results are likely to be more accurate as is seen by comparison 
With an analytical solution. (Table 1). A simpler method to 


locate the spike, reducing the tedious hunting process, has been deve- 
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loped and has been explained in some detail under article 3.6. The 
use of the Runge-Kutta procedure in the entrance is numerically more 


accurate and is also new in elastohydrodynamic studies. 


3.8. Non-Dimensional Parameters 

Dowson and Higginson [59] have formulated three non-dimen- 
Sional parameters, viz. load parameter, speed parameter and materials 
parameter for isothermal elastohydrodynamic contacts. In terms of the 


notation used in this study they are 


Load parameter eee = (330) 
u_U 

Speed parameter SPoq= = 339) 

Materials parameter MP = 4E C3e 10) 


It has also been shown [59] that the effect of load on film thickness 
is very little. For ordinary elastohydrodynamic contacts the variation 
in the materials parameter is also small. Thus, the speed parameter 
emerges as the dominant variable in controlling film thickness. In the 
present study both the load parameter and the materials parameter have 
been kept constant. The effect of speed parameter on the film thick- 
ness has been studied over the range SP varying between 107!2 to 107!°, 
It may however be pointed out that the same numerical scheme can be 
adopted for studying the effect of the other parameters. By varying 
the slip for the same effective rolling velocity temperature effects 


have been studied which form the subject matter of Chapters IV. 
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3.9. Results and Discussion 

The input data and the isothermal results are given in 
appendices E and F respectively. Dowson and Whitaker [19] have pre- 
sented a minimum film thickness diagram showing the three regions en- 
countered in cylinder lubrication. This diagram has been reproduced 
in Fig. F.1 and the range of the present study has been marked therein. 
It is seen that the present investigation is inside the elastic range 
and that the ranges of speed covered are between moderately low and 
moderately high. The pressure distribution and the film profiles have 
been obtained for different rolling velocities and are presented in 
Figures F.2 and F.3 respectively. Referring to Fig. F.2, it is ob- 
served that with the increase in rolling velocity the pressure 'spike' 
moves towards the line of centres. The magnitude of the pressure at 
the 'spike' increases with rolling velocity though not proportionately. 
At low rolling velocities the pressure distribution is closer to the 
Hertzian pressure distribution. 

Referring to Fig. F.3, it is noted that the film profiles 
are parallel in the high pressure region with a depression developing 
in the down-stream end. The area of the depression increases with 
increase in rolling velocity. The film thickness increases appreciably 
with speed but the rate of variation tends to diminish at higher 
speeds. All the above results are in good qualitative agreement with 
Otheserestil tsa (4. lol (ee los? beac 

Cheng [65] has indicated the three thicknesses in common 


use in elastchydrodynamic work. These are the nominal, centre-line 
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and minimum film thickness and are indicated in Pree iy 'Chene [65] 
has also given a table of the isothermal elastohydrodynamic film 
thickness formulae based upon a correlation of the Yresults of various 
authors. A good quantitative check on the results of the present 
study will be to compare the results on the film thicknesses with 
those of other authors, as per Cheng's [65] table. This has been done 
in Fig. F.4 to F.6. Experimental results due to Sibley and Orcutt 
[45] and Crook [47] have been superimposed. It is clear that the 
present studies correlate favourably with experimental results. Fig. 


0.66 
F.6 also yields a correlation of the form hea U 5 


3.10 Conclusions 

The Dowson-Higginson [59] formula or the Grubin [11] formula 
have a narrow range of validity. They yield good results for veloci- 
ties not exceeding 400 inches per second and for Hertzian maximum 
pressure not exceeding 100,000 psi. Also the materials parameter aE 
has to be high in these formulae. This is a severe restriction. At 
low values of the materials parameter the thicknesses predicted are 
much lower than is actually the case. For example both these formulae 
predict zero film thickness for an isoviscous lubricant. This is absurd 
as Herrebrugh's [20] analysis and the measurements of Kannel et al 
[67] have shown. Thus a rapid and accurate numerical scheme becomes 
necessary to calculate the film thickness for all ranges of the 
materials parameter. The numerical scheme devised for the present 


study is good for this purpose. 
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The favourable correlation with experimental work and 
theoretical results also shows that accuracy at every stage is important 
and balancing of the elastic and hydrodynamic film thickness has to be 
done carefully. The closer agreement also shows that the errors in 
the numerical scheme keep within agreeable limits and justify the 


assumptions made in developing the scheme, 
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CHAPTER IV 


THERMAL ELASTOHYDRODYNAMIC PROBLEM 


4.1. General Approach 

The thermal elastohydrodynamic problem requires the simulta- 
neous solution of the Reynolds equation (2.39), the energy equations 
250) (2.52) .and (2.53) along with the elasticity vequations(2 243), 
taking into account the variation of fluid properties given by 
Povattonsa()50)5 and.(2057),.. Asia sclosed.form solution of these 
coupled non-linear integro-differential systems of equations could not 
be found suitable numerical methods had to be used. The method adopted 
here is to use the results obtained in Chapter III as the initial input 
and devise a cycle of operations that will take into account the 


thermal effect on the system. 


4.2. Simplifying Assumptions 

In order to obtain a solution to this highly complicated 
problem, it is necessary to make some simplifying assumptions. Cheng 
and Sternlicht [21] made an assumption of a mean viscosity across the 
film. Cheng [23] removed this restrictive assumption but his form of 
the Reynolds equation has been found to be in error (Chapter II). 
Dowson and Whitaker [24] have used thermal factors, considered functions 
of x to obtain property values at any given X location. This again 
is an assumption of mean Broner vein ieet Further, all the above 


authors [21, 23, 24] have used a surface energy equation as the 
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boundary condition for the energy equation of the lubricant. The 
assumption that is inherent in deriving the surface temperature equa- 
tion is that the problem can be treated as a case of linear heat flow 
[24, 23, 29, 60] and that the dimensionless parameter 


(a) (€ 
Tey 2B 
2%. kK, 





>a Os 


This is generally true in high speed contacts, but in the case of low 
speed, light-weight highly conductive contacts this condition may not 
be satisfied. A less restrictive boundary condition would be to 

equate the heat fluxes at the fluid-solid interface and this has been 


done in the present work. Also, this method yields the temperature 


distribution in the solids. 


4.3. Numerical Method for solving Energy Equations 

A finite difference formulation has been adopted for solving 
the energy equations of the solid-fluid-solid together as a unified 
field. The grid diagram and the finite difference formulae to be used 
are shown in Fig. 6. With the finite difference computational mole- 
cules adopted, the energy equations could be written as a system of 
algebraic equations as shown below. The energy equation (2.52) is 


written as 
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FIGURE 6 Grid Diagram 
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(4.4) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 


(4.9) 


(4.10) 


(4.11) 


(4.12) 
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The boundary condition equation (2.5lc) is expressed as 
K = as 
a = Se = . 
Woipy +25 -2., - Se H]4 aC (4.15) 
= Ay, K , Dy Ke te NEL 
The energy equation (2.53) is expressed as 
= x = 4AX 
ye ptediz eit, (3RS, +) - 
a ape Oe i,pp 2 
o 2Ax = 
T ae FORS eT: = 0 (4.16) 
-1 2 i-2 
= pp-1 ¥5 z »PP i »PP 
Ne PP-1 


From equation (2.4) and (2.5) one obtains the normalised velocity 


gradient and velocity. They are 





(4.17) 
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= _ _wh? 6 Pinan, YM pai. 
yu UaL dx ‘ LU U e chai 
Ee a ee ee 
fe (4.18) 
2u_UaL d& tr <o U 
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where 7 is an auxiliary non-dimensional y-coordinate. Equations 
(e-3) (4.8), (4.9), (4.15) and (4.16) constitute a system of non- 
linear algebraic equations and the following procedure is adopted to 
solve them numerically. 

For the solution of the energy equation the pressure gradients 
and the film thicknesses are considered known [21]. For the first 
cycle of calculations the values calculated from the isothermal solu- 
tions (Chapter III) are adopted. (In the subsequent cycles, the 
pressure gradients are amended by the iterative solution of the modi- 
fied Reynolds equation (2.39) along with the elasticity equation 
(2.43) and the temperature as calculated in the previous cycle are 
used.) Further, for the starting of the iterative solution, the 
density and viscosity of the lubricant are considered as functions of 
pressure only. This reduces the system of energy equations to a tri- 
diagonal linear system and rapid and efficient algorithms are available 
to solve this system. Thomas algorithm [66] was adopted. The adoption 
of this method reduces the storage difficulties and the grid can be 
made as small as is desired. The solution can now proceed colum- 
wise. As soon as the temperatures have been calculated in a particu- 


lar column, the property values are recalculated, adopting the 
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presently calculated temperatures and once again the system is solved 


for that column. The iteration is stopped as soon as the criterion 





on (4.19) 


was met. It was found that in high speed cases acceleration parameters 


[68] of the form 


Go Se (eigen, ee (4.20) 


= 0.47 to 0.61 


(S 
I 


were useful. The convergence was slower in the spike region. The 
Newton-Raphson method was also used in some cases and this also in- 
proved convergence [68]. The block diagram illustrating the numerical 
procedure described is presented in Fig. 7. Simpson's rule was used 
to calculate the integrals in the equations (2.38a to 2.38c) and in 


equation (4.18). 


4.4. Modified inverse hydrodynamic Relation 

It is now necessary to find the changes induced in the 
pressure gradient and in the film profile due to the modifications in 
the temperature field. The methods detailed in chapter III can now 
be used with a modified inverse hydrodynamic relation. Equation (2.39) 


can be rewritten as 
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For discreet points in region B (Fig. 3) the values of ti; Yee t,, 
can be calculated from our knowledge of temperature distribution. 
Further, we can also calculate from the determined temperatures at 
the down-stream end the values of une and Bes Since the pressure 
gradients at both the line of centres and the down-stream end are 
zero, expression (4.21) applied at these points yields 
(- AW.t A) + = (- Mgt oh) | ae (Ce PR) 

where subscript o and * indicate evaluation of the quantities in 
the brackets at the line of centres and the down-stream end respect- 
ively. Equation (4.22) is a first order relation in h* (remembering 
our basic assumption that the centre-line film thickness is insensitive 
to change of temperature) and hi” can be evaluated. The procedure to 
be followed hereafter is identical to that described in articles 3.5 
and 3.6 with the modifications indicated below. 
1. The inverse hydrodynamic relation to be solved is given by 

equation (4.21). 
2. For finding the position of the spike, computer time is saved 

if one starts fromthe location of the spike obtained from the 


previous cycle. Usually only a slight shift towards the line 


of centres is necessary. 


4.5. Further Temperature Calculations 
The temperature field needs to be checked on the basis of 


the modifications in film profile and pressure gradient. One repeats 
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the procedure detailed under article 4.3 with the presently found 
values of the pressure gradient and h. If the presently calculated 
temperatures do not differ by more than 1% of the values obtained in 
earlier cycles, the calculated values are accepted as the solution of 
the problem. Convergence in the temperature field was faster than 

in the case of pressure and film profile calculations. However, 

3 to 5 full cycles were found enough for the range of the present 


work. 


4.6. Calculations for the Frictional Force 
Following Cheng [23], the friction factor is defined as 
the ratio of the frictional force to the load pereunt tewildtisofethe 


roller. This may be written as 


Wo 
ah du ee 
f.== — L d& (4.23) 
rT ul, pee 
h 0.5 ca n 
eid (ae ars 
aes dx 


i 
Hak he adie f = dk (4.24) 


f =—-——L1Jf = dk (4.25) 


Simpson's rule was used for integration purposes. 
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CHAPTER V 


DISCUSSIONS AND CONCLUSIONS 


5.1. Discussion of the Results 

Using a modified Reynolds equation for a lubricant whose 
density and viscosity are functions of pressure and temperature, the 
thermal elastohydrodynamic problem between two rolling and sliding 
cylinders has been investigated. The numerical scheme developed 
analyses the 'isothermal system' and the results serve as an input to 
the thermal problem. One of the contributions of the present work 
is the derivation of the correct thermal form of the Reynolds equation. 
Also, instead of using a surface temperature equation [24, 23, 21], 
the energy equations for the solid-liquid-solid combination has been 
solved as a unified field. This approach extends the validity of the 
solution and has been discussed in Chapter IV. 

A detailed discussion of the numerical scheme for the 
isothermal problem has already been presented in Chapter III. The 
reasons for adopting a tri-diagonal matrix solution for the temperature 
equations and the iterative methods used have been detailed in 
Gneerer EV 

In the present chapter, the results of the 'thermal' 
problem, which are presented in Appendix G, are discussed. The slip 


ratio which is defined as 


Mee UL, 


2 
Slip Ratio = ———— et) 
UL 
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has been varied from zero to a moderately high value of 0.25. 

Figure G.l relates the position of the pressure spike with 
the slip ratio and rolling velocity. Also indicated in the same 
figure is the magnitude of the pressure at the spike. It is observed 
that at constant rolling velocity the position of the spike is moved 
towards the line of centres as the slip ratio is increased. This 
result is in agreement with the result of other workers [21, 24]. 

It is also seen that the magnitude of the pressure at the spike is re- 
duced with the increase of slip. This finding is in agreement with 
the conclusions of Dowson and Whitaker [24] and is in conflict with 
the results of Cheng and Sternlicht [21]. It may be observed that 
with the increase in slip the temperature of the lubricant increases 
with corresponding decrease in viscosity. Herrebrugh [20] did not 
find any spikes for an iso-viscous lubricant and came to the con- 
clusion that spikes are due to increase in the absolute viscosity of 
the lubricant. Thus, any reduction in viscosity should reduce the 
Magnitude of the spike and the present results conform to this 
reasoning. Figures G.2 (a and b) and G.3 (a and b) represent the 
temperature distribution across the film at different x locations 
for rolling velocities 400 in/sec. and 30 in/sec. respectively and 
for a constant slip ratio of 0.25. The presentation of each figure 
has been subdivided into two parts a and b in order to achieve 
clarity. Thus, the figure G.2.a represents the temperature distri- 
butions up to the position of the spike and G.2.b shows the tempera- 


ture distribution from the position of the spike to the down-stream 
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end (x = 0.5). Similarly for the figures G.3.a and G.3.b. Referring 
to figure G.2.a, the temperature rise is seen to be moderate eg Gal 
about xX =- 0.3. Thereafter the temperature increase is more rapid 
and this is attributed to the increase in viscosity in the high 
pressure region. The temperature rise near the fluid-solid interface 
is however moderate. The significant increase in temperature at the 
region of the spike is evident where the non-dimensional temperature 
rise is about 0.185 between ¥ = 0 to x= 0.06. This is due to the 
rapid increase of viscosity due to high pressures at the spike region. 
The mid-stream temperatures are seen to fall beyond the spike (Fig. 
G.2.b). Here, the lubricant temperatures near the walls are seen to 
increase slightly. Both figures, G.2.a and G.2.b, show that the 
viscous heat generated at the centre across the film is much higher 
than that generated near the boundaries. Figures G.3.a and G.3.b 
follow a slightly different pattern. The temperature rise is moderate 
till x= - 0.2. This is attributed in part to the low rolling 
velocity (u = 30). The temperature rise is less rapid in this case 
and at x = 0 the non-dimensional mid-film temperature is 1.14. The 
mid-film temperature then falls slowly to 1.09 at location xX = 0.32 
(not shown in figure G.3.a) due to the fall in pressure. The spike 
occurs further away from the line of centres in the case of low 
rolling velocity. Thus, the effects of decompression cooling and 
conduction are more evident. At the region of the pressure spike the 
temperature rises rapidly. Maximum mid-film temperature occurs at 


the spike. Beyond the spike (Fig. G.3.b) the mid-film temperature 
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Paibissetid) e-soe= 0N5. the temperature across the film is almost 
uniform. A comparison of the figures G.2 and G.3 reveals that the 
solids with a lower rolling velocity reach higher temperatures for 
thetsame slip ratios ~ This finding is in agreement with Cheng and 
Sternlicht [21]. It is also observed that in both cases represented 
by figures G.2 and G.3, the lower solid which moves more Slowly than 
the top solid attains the highest temperature. This result is in 
agreement with other results [24, 21, 23]. Checking the temperature 
distribution at the down-stream end x = 0.5 it is observed that 


the variation of temperature across the film is significant at 


eI 
ul 


400 in/sec. Thus, the assumption of Cheng [23] that t, = 1 at 


i 
i] 


0.5 is inaccurate for high velocities. 

Figures G.4 and G.5 represent the mid-film temperatures 
along the contact for several slip ratios at U = 130 in/sec. and 
U = 200 in/sec. respectively. The temperature profiles follow the 
pressure distribution and maximum mid-film temperature is at the 
position of the pressure spike. In both cases, for zero slip the 
temperature rise is negligible, showing the insignificant effects of 
compressive heating. — For other slip ratios the temperature rise is 
moderate till x = - 0.3 and thereafter the temperature rise is rapid 
till the line of centres is reached. For rolling velocities U = 130 
in/sec. (Fig. G.4) we notice that for low to moderate slip ratios of 
1%, 54 and 10%, there is a slight cooling till the spike region is 
reached. This is attributed to the fall in pressure and hence fall in 


viscosity and the decompression cooling effects. At higher slip 
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ratios of 15% and 25%, the temperature remains constant till the region 
of spike is reached. This is attributed to the fact that at higher 
slips the spike moves towards the line of centres and thus the 

pressure falls slower than in low slip cases. The temperature falls 
after the spike and at low slip ratios the temperature at the down- 
stream end is moderate. As was to be expected, the temperatures are 
higher at the down-stream end for higher slip velocities. A similar 
pattern is followed for the case U = 200 in/sec. (Figs G.5), Hicure 
G.6 shows the maximum rise in mid-film temperatures reached at various 
Slip ratios for different rolling velocities. It is observed that for 
a constant slip ratio the temperature rises steadily with rolling 
velocity till about U = 200 in/sec. Further increase in rolling 
velocity causes the temperature to rise much slower than before. 

Fig. G./ shows the maximum film temperatures plotted against slip 
velocity for different rolling velocities. It is observed that for 

low values of (U, - U,) there is a rapid rise in temperature. Further 
substantial increases in (U, = U,) causes only slow rise in temperature. 
This is attributed to the rapid fall in viscosity at high temperatures. 
Also it is seen that the increase in rolling velocity causes slight 
increase in temperature, showing the major role played by slip velocity 
in the generation of heat. Fig. G.8 shows the temperature rise above 
ambient reached by the lower solid at the position of the spike for 
various slip velocities. At zero (and near about zero) slip velocity 
for U = 400 in/sec. the temperature rise is 4°F and for U = 30 in/sec. 


the temperature rise is negligible. As (U, - U,) is increased, it is 
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observed that the temperature rise is more rapid in the case of lower 
rolling velocities than in higher rolling velocities. This leads to the 
conclusion that thermal stresses are likely to be less severe at 
higher rolling speeds. When it is remembered that the dip in the 
down-stream end is more pronounced at lower rolling speeds, one may 
conclude that if fatigue and thermal stresses are considered important, 
they should be provided for at the lowest rolling velocity the system 
is likely to encounter. 

Fig. G.9 shows the surface temperature rise above the 
ambient for the lower solid along the length of the contact for 
U = 30 in/sec. for different slip ratios. It is seen that the tempera- 
ture rise at the entrance to the high pressure region is of the order 
of 1°F, The film thickness has been found to depend on the viscosity 
of the lubricant at the temperature of the surrounding solids at inlet. 
The insignificant rise in temperature at inlet shows that film thick- 
nesses calculated on the isothermal assumption are valid for the 
thermal case as well. Fig. G.10 shows the plot of friction factor f 
against (U,, - U,)- Super-imposed are the experimental plots of Crook 
[49] and Bell and Kannel reported in the discussion in Reference [23]. 
It is seen that friction factor increases initially with increase in 
slip velocity and then decreases. This is due to the fact that the 
friction factor depends on uy (U, - U,) and initial increase of 
(U, - U,) increases the friction factor. As (U, - U,) is further in- 
creased, To decreases due to rise in temperature and the product 


is (U, - U,) decreases. The trends predicted in the present work and 
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in the experiment are the same, Quantitatively the agreement is 
moderate. The discrepancy may be due to an overestimation of the 
viscosity near the peak pressure. The assumption of a Newtonian 
fluid near the vicinity of the peak tends to overestimate the peak 
as the isothermal analysis of Chow and Saibel [38] has shown. Cheng 
[23] has reported that a slight change of the temperature viscosity 
exponent § causes a significant change in the predicted friction 
factors. 

The problem being non-linear, analytically there is no 
general proof of the uniqueness of the solution [66]. Lf comparable 
results are obtained by using two different numerical procedures, the 
solution is considered correct. This was done in the present study 
for u = 30 in/sec. The film thickness results of the present study 
which compare well with experimental results also show the consistency 
of the numerical solution. 

Assessment of error is very difficult for the problem under 
study. However, conservative error estimates can be made. Pressures 
in the inlet region have been calculated by fourth order Runge-Kutta 
method with a step size Ax = 0.1. Computation of the non-dimensional 
pressures have been carried out with computer errors less than one 
percent. In regions B, C, D the numerical procedures are different 
but the non-dimensional pressures have been calculated with computer 
errors less than one percent. The temperatures are calculated to an 
accuracy of better than ninety-nine percent in all regions. The nun- 
erical procedures used are stable [66] and it is well known that if 


consistency is satisfied, then stability of the system also ensures 
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convergence, 


5.2 Conclusion and Scope for further Work 

The pressure, temperature and film thickness between two 
rolling and sliding cylinders have been determined by solving numeri- 
cally the coupled Reynolds, elasticity and energy equations. The 
isothermal problem was solved first and an iterative procedure was 
developed to solve the thermal problem. The results show that a 
pressure peak exists in the down-stream end of the film for all 
heavily loaded cases. The effect of temperature was to reduce the 
pressure peak slightly in all cases. Increase in rolling velocity 
shifts the location of the pressure peak towards the line of centres. 
Increase in the slip velocity also moves the peak towards the line of 
centres. At low speeds the surface temperature was higher. The slip 
velocity has been shown to be the major factor in the temperature rise. 
Friction factor increases initially with sliding velocity and later 
decreases. The present results compare more closely as regards the 
film thicknesses with experimental findings. The comparison of the 
friction factor with experimental findings shows moderate quantitative 
agreement. 

The present numerical procedure can be effectively used to 
investigate the effects of the material parameter and the load para- 
meter. The effects of non-Newtonian behaviour of the lubricant is 
worthy of study. Study of the dissipation of heat for low reference 
solid numbers Rs, and Rs» will be of interest in thermal elastohydro- 


dynamic studies. 
















mre Tt 


} 
MoW Tian sda ayes One eer) ad ia 
- ‘ e 


rowed s29nipind @i.7 Ore td BY 9e/ma 3 , TWA ont 


1s 
ny pir) f nimyoI9sh 195° ' wysn 2 texts Ya guibite be “aL, ™ 1 
tytooe voyons bia viletaents ,siileaved osiavas ot? d 
; Pike ) 


a 
rT. 7 


tticusti né& bet Sest> beviler ane meldutg Lanse 


oy oft noalaevq Liberendd aid syle oF Sa 


t- - ' 
oF tf mit 4 bis mo“*le-fwoh si? nL avalus 2neq @7aee 


ey . add | fostte eft ,eeeno beheol qaas 


, a 
‘< RLY ! i? vidsaeria dowry Bet 260: 
’ nT ; EWG J Tie Teas 24th a] t aeta nook on3 est 
7 ‘on ueis vViisebses alle of9 &P Saeeee 
) d ao sts er lon 3° ad sisiageh vol FA «4 
’ 7 — 
igJ251 Yolen Sia ot of ovate peed. San See 


bere v2 + wulbtle datw yUletitel geseewenk. yesoa? 


eri t Tite & VISeod FIO STegiTiGS ei Use Jireha tg ant 3938 aa 


162 onl. -<epebeiy.1é im tyes’ haw sscusmin tad 


(i@ @antbai? :utwekit2rsges ‘Watt sedan nots 


“jeu 


:! a v4 








: eeu. wis MwIlts sd o69 S68ungso74 [asta ei ——— - 
~Briiq thes! aie bose iotaneceh tetas bil), 26 stale os #38 -* 


Ly 


4 rs 
BL Aueoeeyd ai ts vod vedsd ndice syalina i aasette ne 7” 












Hos wo ae 202i 40,4 as snares -busi sii edt 
a>? cant aie Pott any as fay fo Se - av tas 
7 | a tbe 7: a en of 2S a é 


. 


ile 


7h 


BIBLIOGRAPHY 


Hertz, H., "Ueber die Beruehrung fester elastischer Koerper", 1881, 
"Ueber die Beruehrung fester elastischer Koerper und ueber die 
Haerte", 1882, reprinted in Gesammelte Werke von Heinrich Hertz, 
Vol. I., pp. 155 - 173, 174 - 196, Leipzig, 1895, English trans- 
lation in: Hertz, H., Miscellaneous Papers, translation by Jones, 
D.E., and Schott, G.A., pp. 146 - 162, 163 - 183, Macmillan, New 
York, N.Y#, 1896. 


Crook, A.W., "Developments in Elastohydrodynamic Lubrication", 


Journal of the Institute of Petroleum, Vol. 49, No. 478, Oct. 1963, 


Pp. 20% 307. 

Martin, H.M., "Lubrication of Gear Teeth", ENGINEERING, Vol. 102, 
Pero. Loo: 

Peppler, W., '"Druckuebertragung an geschmierten zylindrischen 
Gleit- und Waelzflaechen", V.D.I. Forschungsheft, 1938, p. 391. 
Meldahl, A., "Contribution to the Theory of the Lubrication of 


Gears and of the Stressing of the Lubricated Flanks of Gear Teeth", 


Brown Boveri Review, Vol. 28, No. 11, 1941, p. 374. 


Gatcombe, E.K., "Lubrication Characteristics of Involute Spur Gears 
- A Theoretical Investigation", Transactions of the A.S.M.E., Vol. 

OF, 9435 pe Li. 

Hersey, M.D., and Lowdenslager, D.B., "Film Thickness Between Gear 

Teeth, -TransactionscofithevA.S.MsE., Vole */2s 1950 sep. 10a. 


Cameron, A., "Hydrodynamic Theory in Gear Lubrication", Journal of 


the Institute of Petrol, Vol. °38,° 1952, p. 614. 


i rt ; s a. ie : - 
i aR. 
_ - oe 


¢ te 7 4A</ - 2p AG L389 adi” c@ M OW ada 


2015 do kekat "ag Hey tat tos fsstondp A, - 
_ “ Us 
























VittaqooT gare 4 
; @ ‘ 

wy. 
ls ioten) goupfenesh «th sade" | Aes 


7. 
= 

ine ticwlSs q63ee1. eiundapiet off toda 
. Je 


ie vi) Sout wi tedA? 7 , eet , sree - : 


Ley had ACs aya pert - €€. og gel ghee 
jot puter focakt Of ate cal eanl 

f —@SD sie ,.A.o <stones i a 4 
elas te treY 

: thi a level ,.WcA ,dleed as 

ae 


WEboLIS? 20 .63/ Taeal sis Sp Lasapot 
P ve. = 205 ‘Tt 


y 07% 4g ad 

7 iu Sie’ Iecerdouere ng erabqaat 
[Saree 4509.7 "ne ose lzalegW fey —sield 
HAT wil os selsuseee ne! Gee litebteat of 
éy Sua arr 34 ai 20938 sc? to bow e188 


Ol £0) a4) 62 .i6¥ wise Saved sya. 
i tison ga) 4100 tase aizdut"’ ae Pe is tua 


= 


ON ag tk 48 7 
eet ee a 


ice 


— = 





Poe 


as 





Su t aA « “ 


7 





72 


9. Mcewen, E., "The Effect of Variation of Viscosity with Pressure on 
the Load Carrying Capacity of Oil Film between Gear Teeth", Journal 
ofptheminstitute ofsPetrol,sVolon38. 01952) 39646. 

10. Blok, H., "Discussion. Gear Lubrication Symposium. Part I, The 
LubricationgoteGears!s J.ilnsteiPetrol savolsn38eebas673% 

11. Grubin, A.N., and Vinogradova, I.E., "Central Scientific Research 
Institute for Technology and Mechanical Engineering, Book No. 30", 
Moscow, 1949, DeS. IT. ReeTranslation Nowe337; 

12. Weber, C., and Saalfeld, K., "Schmierfilm bei Walzen mit Verformung", 


Zeits. ang. Math. Mech., Vol. 34, 1954, p. 54. 


13. Dowson, D., and Higginson, G.R., "A numerical Solution to the 
Elastohydrodynamic Problem", Journal of Mechanical Engineering 
Science, Vol cas ENe. 1s 1959. fppeaGe-215. 

14. Dowson, D., and Higginson, G.R., ''The Effect of Material Properties 
on the Lubrication of Elastic Rollers", Journal of Mechanical 
Engineering Science, Vol. 2, No. 3, 1960, pp. 188 - 194. 

15. Archard, G.D., Gair, F.C., and Hirst, W., "The Elastohydrodynamic 
Lubrication of Rollers", Proc. Roy. Soc. Series A., Vol. 262, 1961, 
ppe,5it=an 72% 

16. Archard, J.F., and Kirk, M.T., "Lubrication at Point Contacts", 
Proc. Roy. Soc. Series A., Vol. 262, 19GIRepp. soe =O Or 

17. Dowson, D., Higginson, G.R., and Whitaker, A.V.,"'Elastohydrodynamic 
Lubrication: A Survey of Isothermal Solutions", Journal of 
Mechanical Engineering Science, Vol. 4, No. Ziwl962¢e pp .W121gen 126. 

18. Rhoads Stephenson, R., and Fletcher Osterle, J., "A Direct Solution 


of the Elastohydrodynamic Lubrication Problem", A.S.L.E. Trans- 



















Fania w apiiiisel % goeria. aft? eed oe anh 


| 7) 
ea Dy asavuiad mpi Fie V2 S28) |aneviT eo paod fs ae 


7 

“o0L 60 fev, hae ae Atul esoms ails 40 ’ 

' a aa 7 

a axe . vtisotvénd saga (aeltagoel” _ ib paler 
i 


‘Ja7' (gah) ot , eta Fo soliasipG@al 


3 i ehal #yvoberaceiVv Oae...%.d pai 
hain Ber ‘epo lori! oy 74) 930395 


biplongut-iais Asa «Shes woopait 

is bunfaee” A By i Luc Gar , od yredaw 
é Cty POT po iaatl pte mis AOR ok | 
; rau A’ 66a) | doeeigett base .<0 ,penweie 


él 
, 
Edo} . Vay? bot yaigg hake 


7 — — i ee 7 


J € . oa ts 


7 | ee i; P HUBBES os) ae ’ »4] ‘ roawad af it : 
; . i 4 7 


¥ [at «\roeeta: te geben tale ald ne n 


ir el iS «Let eae vii sap hee 
ea ' ise 


ov? Soe, SLT See Oe Janik - 


i 
}a4 Mod iB (oclug #9 maeretey 4] 
: 
ee a | pe 
| iA8° mfsentyddl” VATA G seh “tem tt aan a 
A) ae ae ee ous vat T. ae 
, Serarereimny 24h la”, A pee — fins” 4M, eesaeagllt ig A *¢eaiawatt K i 
ol, aacrryL OR Temgagsned 26° way a 
x a oa | ma 


| on 







Bh ist Le 
a 7 


i 
cat i , 2 4 













WY ., aa 


she 


20. 


21. 


Fane 


23. 


24. 


25% 


26. 


(te) 


actions, Vol. 5, 1962, pp. 365 - 374. 

Dowson, D., and Whitaker, A.V., "The Isothermal Lubrication of 
Cylinders", A.S.L.E. Transactions, Vol. SIBLIGSS pp. 2224858234" 
Herrebrughyert , "Solving the Incompressible and Isothermal Problem 
in Elastohydrodynamic Lubrication through an Integral Equation", 
Journal of Lubrication Technology, Transactions of the A.S.M.E., 
Jan. 68, pp. 262 - 270, 

Cheng, H.S., and Sternlicht, B., "A Numerical Solution for the 
Pressure, Temperature and Film Thickness Between Two Infinitely 
Long, Lubricated Rolling and Sliding Cylinders, Under Heavy Loads", 
Journal of Basic Engineering, Transactions of the AsoeMcha, 

sept. 1965, pp. 695 = 707% 

Dowson, D., Discussion to Reference 21, Journal of Basic Engineering, 
Transaction of the A.S.M.E., Sept. 1965, pp. 705 - 706. 

Cheng, H.S., "A Refined Solution to the Thermal-Elastohydrodynamic 
Lubrication of Rolling and Sliding Cylinders", A.S.L.E. Trans- 
attions,+ Vol.o8, 019655- pps 3979=" 410. 

Dowson, D., and Whitaker, A.V., "A Numerical Procedure for the 
Solution of the Elastohydrodynamic Problem of Rolling and Sliding 
Contacts Lubricated by a Newtonian Fluid", Proceedings of the 
Institution of Mechanical Engineers, London, Vol. 180, Part 3B, 
Paper 4, 1965, pp. 57 - 71. 

O'Donoghue, J.P., and Cameron, A., "Friction and Temperature in 
Rolling Sliding Contacts", A.S.L.E. Transactions, Vol. 9, 1966, 
pp. 186 - 194. 


O'Donoghue, J.P., and Cameron, A., "Temperature at Scuffing", 















—s Oa0Tf 4 f sar s 
, cA ay ¢ (AOV @® 
V j 4 ¢ i 
t qa AIL ' 4 
au oo 
ii J aed 7 ee awh e < Tener, tly ' _ 


et PL | Be By ht 1 
moans 4Ad ealo tee” cdl dgurtdatealt 0% 


: ‘ -_ 
ve , Ths i Da ey ‘i (od aha ak be 





vt owt ' aotaeattdeld je Isnyuet 
OSS = 40 .9e .8O sane 


=e 1) ' 
srlersye fae , 2.8 ,aaed® Ie 
> , | 
o 1! ; 
[4 SAB reuaeT ,oTeuagess 


toy yhvihs as “Se ae ; 
VER (a Yee We ) yas i... - yet’ ES 

; . bf coti tank 2 . Hsenkadad d 
Lit = V0r RAM (Be hav cetiohton 


HA" AA oo ie 2 ost 
nenyborbvidedaadd add to cots 7 





bs tear ediud 






ath) 


Saety Late bis “ol soung A shore) bers eit pom ay 





Zh. 


28. 


712 


30. 


el. 


Bere 


25% 


74 


Proc. Instn. Mech. Engrs., London, Vol. 180, Part 3B, 1965-66, 


ppwtS5. —B94: 

Blok, H., "Theoretical Study of Temperature Rise at Surfaces of 
Actual Contact Under Oiliness Lubricating Conditions", General 
Discussion, Lubrication, Institution of Mech. Engineers#OVol. «2, 
LOSTSGpp m222 = 235. 

Carslaw, H.S., and Jaeger, J.C., "Conduction of Heat in Solids", 
Oxford University Press, Second Edition, 1959. 

Jaeger, J.C., "Moving Sources of Heat and the Temperature at 
Sliding Contacts", Journal and Proc. Roy. Soc. New South Wales, 
mol. 5265019435 ppee203 — 224. 

Allen, D.N. de G., "A Suggested Approach to Finite Difference 
Representation of Differential Equations with an Application to 
Determine Temperature Distribution Near a sliding. Contact", 
Quarterly Journal of Mechanics and Applied Mathematics, Vol. 15, 
19 62— ppallli=s33: 

Symm, G.T., Cameron, A., and Gordon, A.N., "Contact Temperatures 
on Sliding/Rolling Surfaces, Proc. Roy. Soc. Series A.,Vol. 286, 
pp és45 i=a6 1 

Francis, H.A., "Interfacial Temperature Distribution Within a 
Sliding Hertzian Contact", A.S.L.E. Transactions, Vol. 14, 1971, 
pp. 41 - 54. 

Cheng, H.S., and Orcutt, F.K., "A Correlation between the 
Theoretical and Experimental Results on the Elastohydrodynamic 
Lubrication of Rolling and Sliding Contacts", Paper No. 13% 


Proc. Inst. Mech. Engrs., London, Vol. 180, 1965-66, p.158. 


{ * 

' ® Phe ‘ 

' «) c st hh spelt Rheem Cm 
: rw j 

} i i ne Ul i 


' 
“~~ 
, {4 I 
1 
r 
; - i 
—— 
- at ‘ ’ 
4 ~* €% # 
— 
i : y t* I C ritve na PL 
+ ani j 
fun - —- 
~ al * 
: . staf 
i + 
~v a —te Py 
if + he >> we 







o“ 


iva ineaod 


ie , -2.H) ogued s€ 
+ pad fa 
- a 


od neawsed febsaketind AY , 4 2 sua%0 hb 
— : ' 7 @ i . 2 : 4 
7 


wien ie ae 








75 


34, Milne, A.A., "A Theory of Hydrodynamic Lubrication for a Maxwell 


Baquid'., Proceedings of the Conference of Lubrication and Wear, 
Instn. of Mech. Engrs., London, 1957, Paper 41. 


35. Tanner, R.I;; "Full Film Lubrication Theory for Maxwell Fluid", 
intemational Journal of Mechanical Sciences, Vol. I, 1960, 

Bp- 82064-5215. 

50, Hell #91 ..Fos "Elastohydrodynamic Effects in Lubrication”, 
M.Sc., Thesis, 1961, University of Manchester. 

37. Bell, J.C., "Lubrication of Rolling Surfaces by a Ree-Eyring Fluid", 
Transactions of the American Society of Lubrication Engrs. , 

VOL Wt) sl 962 eps 160. 

38. Chow, T.S., and Saibel, T.S., "The Elastohydrodynamic Problem with 
a Viscoelastic Fluid", Transactions of the A.S.M.E., Paper No. 70 
=sLubJd. 

39. Merritt, H.E., "Worm Gear Performance", Proc. Instn. Mech. Engrs. , 
Wondon, aVol< 1129, 1935e.n.0127. 

40. Lane, T.B., and Hughes, J.R., "A Study of the Oil Film Formation 
in Gears by Electrical Resistance Measurements", British Journal 
of Applied Physics, Vol. 3, 1952, pp. 315 - 318. 

41. Crook, A.W., "Simulated Gear Tooth Contacts", Proc. Instn. Mech. 
Engrs. London, Vol. 171, 1957, pp. 187 - 214. 

42. El-Sisi, S.I., and Shawki, G.S.A., "Measurement of Oil Film Thick- 
ness between Discs by Electrical Conductivity", Trans. A.S.M.E. 
Journal of Basic Engineering, Vol. 82 D, 1960, pp. 12 - 18. 

43. Leach, E.F., and Kelley, B.W., "Temperature - The Key to Lub~vicent 


Capacity”, A.S.L.sE. Transactions, Vol).8y 1965; pp. 271 = 285. 



















hate 4. a0) Tweed oebeworbel Gy Bleert ON" (5A AR =fy4 df 


ee 


7 | 1 20t989 fd Re Jor coma 2 : 
WV. sce 198 hs .teeadl git Ro- Reumnaot. inno) seepage, 
C19 =) ae ae 


v ina ie 
a bit +79 shewieberpyreaen ti”, .Weteg hie 


plieteytwd.,leeL  ebeedl 2 .a8 . 

0 -95¢2)n2 potdte?, To! setwaottdd” pul 

An 

to: vant dO? ties tems ang Wy sie Lae | 

ASL «@ (S00 (Ei vie 

ala iT ee « Ladiee iim , 2868 , wad. 

slat p(T Pith is (GE De ried “bin Ll? os Seeienel? ay , 
a ~*~ . i Lg 

- sidutl coal | 

Wie reTqu teed mov". 21h eae 

sl oy -@fGl ,2tl fet ,webeed 

"i? A” Piva ty angen Stim ‘ .&7T sont 

“A ovovTe tool! Teo heaei>- i) ataietege ae 

me’ .£ OV eptee tt seen 


port's ‘top 19 puesl Tae) bets Sie? oie 


HS ~ Th aq 9 TGRE TUE ot ation acme 


if 
ia 


bE L445 126 3o soeonweeelt (ABO conte hee gi ne 


a Sey izawbaeD inp rr2oeta td moet am 

















3 


= a 
To 7 
hp aE, 7 ' Pas? ’ 


ae : : ie - 3 a se ; 


44, 


45. 


46. 


47. 


48. 


49. 


Sit 


SEs 


52. 


76 


Takitan@ 7 cE.) Chiu; Y.P., Huttenlocher, D.F., Kamenshine, J.A., 
Sibley @1nB.eand Sindlinger, N.E., "Lubricant Films in Rolling 
Contact of Rough Surfaces", A.S.L.E. Transactions, Vol. 7, 1964, 
pp. 109. — 126, 

Sibley, L.B., and Orcutt, F.K., "Elastohydrodynamic Lubrication 

of Rolling Contact Surfaces", A.S.L.E. Transactions, Vol. 4196. , 
pp. 234 - 249, 

Crook, A.W., "The Lubrication of Rollers", Philosophical Trans- 
actions, Series A., Vol. 250, 1958, pp. 387 - 409. 

Crook, A.W., "The Lubrication of Rollers. II. Film Thickness with 
Relation to Viscosity and Speed", Phil. Trans. Series A., Vol. 254, 
196150 pp. 223.- 236. 

Crook, A.W., "The Lubrication of Rollers. III. A Theoretical 
Discussion of Friction and Temperatures in the Oil Film", Phil. 
Erans .sceriess Aye Vol? «25450 pp... 237e— 258. 

Crook, A.W., "The Lubrication of Rollers. IV. Measurements of 
Friction) andsEffectivesViscosity> “Phiivetrans. Seriés@Al ,“Vol.e255, 
1963 ,Tpp.t28k = 312% 

Dyson, A., and Wilson, A.R., "Film Thicknesses in Elastohydro- 
dynamic Lubrication at High Slide/Roll Ratios", Proc. Instn. Mech. 
Enperse, V0Ol.e103."Part OP,sPapers 11, 1968 - 69, pp. 81 = 97. 
Kannel, J.W., Bell, J.C., and Allen, C.M., "Methods for Determining 
Pressure Distributions in Lubricated Rolling Contact", 

A.S .L.E.s Transactions, Vol. 8,°1965, pps.250°= 270, 

Kannel, J.W., "Measurements of Pressures in Rolling Contact", 


Proc. Instn. Mech. Engrs., Vol. 180, Part 3B, Paper 11, 1965-66, 
pp. 135 — 142; 







" 





















tT] 
‘ 7 
oa leaner c Fall sls ale tae 4 | ur e My asthiet 7 - 


Ul 


fina r Serpe ; Te iT.  naey ‘knee Bm, aT elle 
AOL 2% ‘V7 nilisqgunss) .5)t LAy" seo eet giet, to Baa 
Act — Sg 
; > baits siteels r oe oss “— baw ,. 0.5 yeetere 20 
iulsta hy.) @2ob TRat> 3 aM? qutttes ta” 
0st - $6 ae 
tr aed teoeail ye” MA _, 4009) 43h 
Fi cy . SOL’ (Ore sie +f wits , enor? 5a % 
th itd r rect aft » SA ee it 
Pape Lai,» 4 ied bee yiretcalV ef ebfaeiet 
SEL ~ ELE .qq, gee 
tu Vy, ae LIpa ese oT”. Wie eae 


rue not Sve «0oio be Yo cobs 


Leo ~ he y set tay rr | Hoey et 5 ceemre 


‘taille te inlosetvelny oer” , sa _ lees? . ; 


> 


' PLA Ne Pe efV su liaethy) how aolsalet | 
: a 

Sif 18% -og) 2B 

ne" A ie et 


- 


31 [iea\ebil? «eld do ob teeta dis ns 
 - 18 .99 89 ~ 08! Thelen 


Detniaepefn]-<-3? abowsas"” ,.Mo gor n ie sede eek <oWel ‘fata h 









) Sietons gar ti ok bedastedus it anettn@ie ras? some jn 
an oF : = = 5 5 ‘ 





yh 4 » = 


- ~ - 
ae 





ia 


23 


54. 


By ke 


56 


Dil. 


58. 


a9 


60. 


61. 


77 


Longfeld, M.D., "Pressure Distributions in a Highly Loaded 
Lubricated Contact", Proc. Instn. Mech. Engrs. , Vol. 180, 
Party3D. pp. g) 13h. 118. 

Orcutt, F.K., “Experimental Study of Elastohydrodynamic Lubrication", 
A.S.L.E. Transactions, Vol. 8, 1965, pp. 381 - 396. 

Cameron, A., and Gohar R., "Theoretical and Experimental Studies 

of the Oil Film in Lubricated Point Contact", Proc. Roy. Soc. A., 
Vol. 291, 1966, pp. 520 - 536. 

Westlake, F.J., and Cameron, A., "Optical Elastohydrodynamic 

Fluid Testing", A.S.L.E./A.S.M.E. Lubrication Conference, 
Pittsburgh, Pennsylvania, October 1971, pp. 101 - 111. 

Sanborn, D.M., and Winer, W.0., "Fluid Rheological Effects in 
Sliding Elastohydrodynamic Point Contacts with Transient Loading: 

1 - Film Thickness", Transactions of the A.S.M.E., Journal of 
bubrication Technology, April 71, pp. 262 — 271. 

Sanborn, D.M., and Winer, W.0., "Fluid Rheological Effects in 
Sliding Elastohydrodynamic Point Contacts with Transient Loading: 
2a lraction. « bransactions.of the AwS.M.E., Journal of Lubrication 
Technology, July 71, pp. 342 - 348. 


Dowson, D., and Higginson, G.R., "Elastohydrodynamic Lubriaction", 


Pergamon Press, First Ed., 1966. 
Cameron, A., "Principles of Lubrication", Longmans Green and Co. 
Ltd. ,.Pirst Eds, 1966. 


Timoshenko and Goodier, "Theory of Elasticity", Mcgraw-Hill Book 


Coe, Inc. oecOnds nae ton. 








—“ 


' 


7 
°7y ae 


¥ 7 
(i otigesso” . a." ,bheRigee 
" 
; +4 " isa ww batsoltd 

a) 


° 
- *§ 
] 
7 ‘ 
1 ~ _ _ 
bs J0S9 aay 4 a» 
’ ‘ 
~— a 
a ~ an f 
4 
¥ 
im 
' + 
ri 
@ 4 
‘ 
§ 
La - . * 
3 % 
é - a ' al an 


‘ poe «@ A ‘ —_ 
7 wiv ine ,.t.4 ,asedase: 2ae 
cw - 
- ia % b/ fesaets. ealbi le 
af ~ 
; “Sembisid? wld =f) 
< 
gGtopniloa? op2 ipobaded - 
=i 
i3nld. bas, .4.0 jesednee® . 
» @ 
’ rbydodeata gntbele 
¢ f +6 whiobonest -"eottoarT + 8 


“fn f i iA we pay it ° ated —_ iu — o bd as 
! PIBEl? e. to) q 4 teriiget ‘ Vie t «Ol -%O we] whe 

f 3 \ = 
+t Steer «ieee ' : 
‘OO0L , sb 2070T ,eaeeS ome 


o 





oo eT: Tt 4%) =] ’ > ' u : , : 
a ie (fie as ; ave be ' ny aa abe sf" 
ud h rt 4 bi “s . a 6 Pian Bie ay ce: rt 40) SeLQglarr ’ e A Py 





62. 


63. 


64, 


65. 


66. 


Gy 


68. 


78 


Conte, S.D., "Elementary Numerical Analysis", Mcgraw-Hill Book 
Co., 1965, First Ed, 

Wernick ye. J.. Some Computer Results in the Direct Iteration 
Solution of the Elastohydrodynamic Equations", 62-TR38, 
Mechanical Technology Incorporated, Rport, Feb. 1963. 

Poritsky, H., “Lubrication of Gear Teeth including the Effect of 
Elastic Desplacement", First National Symposium (Fundamentals of 
Friction and Lubrication in Engineering), American Society of 
Lubrication Engineers, 1952. 

Cheng, H.S., "Isothermal Elastohydrodynamic Theory for-the full 
range of Pressure-Viscosity Coefficient", Journal of Lubrication 
Technology, Transactions of the Se mvinio ee: Welles CWP soy, ey 7Acy 
Lee, Stanley E., "Quasilinearisation and Invariant Imbedding", 
Academic Press, New York, First Edo tcl ROG. 

Penne lie J. Woet al A Study of the Influence of Lubricants on 
High-Speed Rolling Contact Bearing Performance", ASD Tech. Rep. 
No. ASD-TDR-61-643, 1969, Part VIII, June 1968. 

Hahn, E.J. and Kettleborough, C.F., "Solution for the Pressure 
and Temperature in an Infinite Slider Bearing of an Arbitrary 


Profile", Transattions of A.S.M.E., Journal of Lubrication 


technology. Oct. 1967,) pp. 445 = 452, 


' 





j 
as 
4 


a" |). 2.0 gaedd. Se 


—— a 





= 
‘ - , ‘ 

aol TMP £40 srpleosetast :,' ea lc 

= = 

2 


Cok --qy S0Cl «cad -geeebagisee 







enna t 


hay 


Al 


APPENDIX A 
RUNGE-KUTTA 4TH ORDER INTEGRATION TECHNIQUE 


Consider an equation of the form y' = f (x,y) where the 
prime denotes the derivative dy/dx. The value of y at station 
X41 in terms of its values at x is obtained by Taylor's series 


expansion as 


2 =) 


> h h 
A ahah gh yah hy! + = u as Speen + -------- CARL) 
x FL x +h where h is the step size 


The Runge-Kutta method proceeds from station x to 
station X41 using values of the first derivative only calculated 
at a number of intermediate points. They are chosen to give agreement 


with the first few terms of the Taylor's series shown in (A.1) above. 


A commonly used set is 


Ky = £ (x ya) 
Kjh 
= =f (x, + h/2, Ne epee 
Kh 
es (x, “tony 2 Peas =m 
K, = f (x, +h, va + Kh) 
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The local truncation error of algorithm (A.2) is 0 (h°). The method 
requires the evaluation of the first derivative four times per each 
step but this is not much of a problem in a computer. It has the 
important advantage that it is self-starting, i.e. it requires only 


Eiesvaluesof sy “at a-point =x = x inorder, to find “yo sand jy") at 


In the present study, equation (2.40) along with equation 
(2.41) can be solved directly through this method if h can be calcu- 
lated. Im the first cycle of calculation, h can be calculated from 
the Hertzian relations (2.27 to 2.30a). In the second sycle of calcu- 
lations h is calculated at discreet points through equation (2.26) 
and the intermediate values of h needed for the Runge-Kutta technique 


is obtained through an Aitken-Lagrange interpolation. 
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APPENDIX B 
CALCULATION OF ELASTIC DISPLACEMENT 


In the calculation of the elastic displacement of the 


cylinders we encounter an integral of the form 
B 
I (x) = fp (&) Inlé-x|de (B.1) 
-A 


If care is not exercised in the evaluation of (B.1), the singularity 
in the integrand can cause large errors. The singularity is removable 
by adopting the technique of approximating the function p by a poly- 
nomial in each sub-interval, performing the integration in closed form 
in the sub-interval and summing over the whole region of integration. 

The choice of the degree of the polynomial approximation to 
p is arbitrary depending on the nature of p. In the literature both 
linear approximation [18] and zero degree approximation [15] have been 
used. However, since the pressure function can have considerable 
curvature, the method developed by Wernick [63] is adopted. 

Wernick [63] subdivided the interval -A to B into N 
sub-intervals requiring only that the value of p (&) be available at 


the end points and mid-point of each sub-interval. He represented 


I, (x), (B.2) 
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where 


on 

1, (x) = [4 p (ein|e-xlac (B. 3) 
ap 

J 


and obtained 


I, (x) = [(E-x) {In| &-x|-1} {p(e)—-5(e-x)p'(E) + 


(E-x)* p"(E)/6} + (E-x)?2 {p'(E)-5(E=x) p" eye 
p (2)=5(E—x)p (2) /93/41 (B.4) 


5 
where 
' =p! eae Ope 4 De -?D, Zi. 
ae (E,) ( P, Pia, Py4)/ A, 
Error 
' 1 = a = 2 
epee? Gay (P, aes oe egg ) — ) (B.5) 
: ) 
it Ww ) 
Boer Be (P, 2 Pas F P54) /A2, 
and 
Kees : - &, Be6 
Be ecg ta (B.6) 


This method can be further extended to yield the pressure distri- 
bution if the deformation is available and the pressures at the two end 
points of the interval where the deformation is available is known. 

A series of linear simultaneous equations can be set up and solved for 


the pressure distribution. This has been called in the present study 
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as the ‘inverted Wernick’ method. It will be noticed that the sub- 


intervals need not be uniform. A comparison of the results obtained 


through Wernick’ 


in Table of. 


Ss method and an analytical solution has been presented 
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APPENDIX C 
ITERATIVE METHOD FOR CALCULATING POLYNOMIAL ROOTS 


Given a polynomial 


£2) = 9s az (C.1) 


Let z= x+iy be a starting value for a root of f(z). Then 


z= (x + iy)” (C.2} 


Define x, as real terms of expanded equation (C.2). Define y, 3s 


imaginary terms of expanded equation (C.2). Then for 


n= 0 
Xo = 1.0 
es 0.0 
n> 0 
a Oe ye CGe3) 
Yn ~ *Yn-1 a y*n= 1 (C. 4) 


Let Ube the real terms of (C.1), V be imaginary terms of (C.1) 
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N 
Usaaeey ax, (€G..5) 
n=0 
N 
Vea Daan (C.6) 
n=0 
or 
N 
U = ay + 2 aX (Coy) 
n=1 
N 
V= 3} ayn (G58) 
n=1 
N 
dU 
n=1 
N 
oU 
oy Hache ae eer 


Equations (C.3), (C.4), (C.7), (C.8), (C.9) and (C.10) can be performed 


iteratively for n=1 to n=N by saving X-1 and y Using 


n-l° 


Newton-Raphson method for computing Ax and Ay the result is 
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after applying Cauchy-Rieman equations. Thus, for the next iteration 


x" x + Ax 


“i 
It 


yoo ny, 


This procedure was found accurate for solving the cubic inverse hydro- 


dynamic relation (3.3). 
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APPENDIX D 
CENTRE-LINE SHIFT TECHNIQUE 


The elastic displacement (equation 2.26) can be written 


in the form 


h = hy + fo) ent (Gea) 


Dowson and Higginson [59] have shown that any linear correction to 
the displacement can be done arbitrarily. If the centre-line of 
the cylinders be moved a distance Ac relative to the pressure distri-- 


bution, the film thickness becomes 


= Z 
h' = h + Getey s (D. 2) 
The change in film thickness is ao + constant. Therefore a 


centre-line shift of Ac produces a linear change in the film thick- 


ness, or in other words, a change in the slope of the datum of - Ac/R. 
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APPENDIX E 


INPUT DATA 


Diameter of rollers (Ry; R,) 
Material of rollers 

Load (W) 

Pressure viscosity exponent (a) 


Temperature viscosity exponent (8) 


Pressure temperature viscosity exponent (vy) 


Inlet viscosity Cu.) 

Modulus of elasticity (E,; ES. 
Poisson's ratio (0), o5) 

Inlet temperature (T.) 

Thermal conductivity of oil (K) 
Specific heat of oil sey 

Density of oil (p) 

Thermal conductivity of solids (Ky; Ky) 
Density of solids (o 1 > e,) 


Specific heat of solids (C_ ,C_) 
it | 1b) 


Pressure density coefficient (Cy) 
Pressure density coefficient (C,) 


Coefficient of thermal expansivity (ce) 


3 inches 

steel 

OC32e cons / ins 

- 0.214 in.2/ton 
8300° R 

350) in, 2°F/ ton 
0.4 poise 

13500 tons/sq.in. 
Oras, 

VOTE 

Op leB fee Fein. 
C.4 B/1b.deg.F 
Q203250ibe/in. 
26.69B/i twat. 
Qe 283etbey) te 


On) B/ ibe kr 


0.009 sq.in./ton 
0.026 sqcin. /ton 


0.000348 per °F 
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APPENDIX F 


RESULTS OF ISOTHERMAL SOLUTION 
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FIGURE F.1 


Intermediate Elastic 





Minimum Film Thickness Contours 
and Range of Present Study 


F2 





F3 


SUOLINGLUZSLG aunssaud 


D/Xx 


00c= 


28S/Ul OOP =/ 


¢°d wqanolA 


Sq| OV9 


M 





NM/oa 


‘$a ONS = we 





| 


F4 


9 


SoTtfoid Witty 


D/X 
cl 80 v0 0 


Of 


(2esut) n 


€°d wanold 


y0- 8°0- cl- 


Sq| OV9=M 





°u/4=H 





F5 


5x104 





10°5 


——-——-=- ARCHARD, GAIR AND HIRST (15) 
DOWSON AND HIGGINSON (17) 

—-—- CHENG AND STERN-LICHT (21) 
PRESENT WORK 


A SIBLEY AND ORCUTT (45) 
@ CROOK (46) 





10°7 owe 10°75 
2p (a+ F-)U ( ee 
R 


FIGURE F.4 Comparison of Theoretical and 
Experimental Minimum Film Thicknesses 
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APPENDIX G 


CURVES AND RESULTS FOR THERMAL SOLUTION 
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Mid-film Temperature Rise for U = 130 in./sec. 


and different Slip Ratios 


G7 





sa2\ni OETEU 
at O he Ww 


0 
-0.6 


FIGURE G.5 
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Mid-film Temperature Rise for U = 200 in./sec. 


and different Slip Ratios 
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APPENDIX H 


FORTRAN PROGRAM 


RESULTS OF THE ISOTHERMAL SOLUTION WERE STORED 


IN LINE FILES AND USED FOR THERMAL SOLUTION 
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THERMAL SOLUTION OF ELASTO-HYDRODYNAMIC PROBLEM H 
* Y MATRIX METHODS 


IMPLICIT RFAL*38{( A-H,0-Z) 

DIMENSION DUY( 116101) .0SX(11,101)s0SY(11,101)eS{(11 
KsLORIS TC LLsSLTON) »sVULLs LOL) sZ1¢(11st1)sGXS8(101) FxXB(1 
*¥O1) ,HC101),HD(101)9,PD(101),PT(101),Z2¢1511)5AC101) 
*ePC101)200101) ,UC115101)5TLC1 (200) 6 TLC2( 200) »sDMUE ( 
* 200) 

COMMON RR,» AL»HOsVSsPSsPMAX.2UD 


INITIAL VALUES ARE ASSIGNED 


READ(55101)U23U1.0SeVS 
FORMAT(4F10e6) 


ASSIGN INITIAL VALUES GF PRESSURE 


NONOIMENSTONALISE 


UD=VU14U2 

CALL HERTZ(HD,»PD) 

CALi THERMACHD.PD,TLCIL,TLC2 SY) 
SToP 

END 
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rand tim y 
7 ev avrg 
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ni 


SUBROUTINE HERTZ(HD»PD) 

IMPLICIT REAL*8(A-H,O0-Z) 

DIMENSION HO(101),HNO(101)5PD(101) »PND(101) »DPXND( 
*(101),00101)0V(101).XCOF(4).,ROOTR( 3) ,ROOTI( 3) sHTTK 
*101)sCOF (4) sDELPD(11.1).PDR(101 )sPSPI(191)e¢HSPI(19 
*¥1)sPDIPCIDLISHTT1(101) sHSPIIC191),HT( 101) 

COMMON Re AsHOsVSsPS»PMAX2U 

EXTERNAL FUN 

HO=0e1620E-04 

READ(5410)PL3sS51sS52.E1sE2sR1R2 

T=(1 e--CS1¥*S1))/E14+(1 e-(S2%*52))/E2 

PI=3214160D00 

EC=T/SPI 

R=RL*R2/(RIFR]) 

A=2 6 *¥DSQRTC(R¥PL#¥EC) 

PMAX=2e*PL/Y(PI*A) 


CALCULATE HERTZIAN PRESSURE AND DISPLACEMENT 

PS=2.25D00 

WRITE(6 s20)PL oe SIL sS2sEl sE2sR1eaR2aRs AHO 
FORMAT(10F13.6) 

WRITE(6525) 

FORMAT(T10.,*HERTZIAN PRESSURE AND DISPLACEMENT®//T 
* 15.*NODE* »T402*P 

ARESSURE IN TONS/SG@QeINe* sT759* DISPLACEMENT IN INCHE 
* S*) 


39 


{ad 
ul 


40 


X=-5Se 
90 30 f1=1,41 
Z1I=(A*¥A)/(2-*R) 
Z2=DSQRT (DABS( X* X-1.2)) 
HDCI 9=Z1 *(DABS(XK¥*Z2)-DOLGG( DABS{(DABS( X)422) ))+HD 
POC 1I)=O0- 
X=X+0.1 
X=—-1le 
9G 35 1=41,61 
POCIJ=PMAX*DSORT (DABS(12-X*X)) 
HD{ 1 )=HN 
K=X4+0e0l 
J=41 
DO 40 1=61,101 
HD(T)=HD( J) 
PDN(1=0- 
J=J-1 
90 50 [=1,101 
WRITE (6s45)1T.PD(1),HD(T) 
FORMAT OTIS es 1L2sT454Fl5e85T75eF1508) 
CONTINUE 


CALL FCAL(ReEC HOLA SPD HT 1019101 299.,HD.1) 
CALCULATE OPTIMUM FILM THICKNESS 

A RUNGE-KUTTA ALGORITHM IS USED 

ITER=1 
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70 


64 


80 


61 


90 


120 
110 


190 


ste Te | 


130 


H 4 


CALL RRK2CFUNs 06 12-50 0900001940sPND) 
DO 60 1=1,.41 

POUT )=PNDC I) *PS 

CONTINUE 


CHECK WHETHER ASSUMED VALUES OF MINe THICKNESS IS 
x OeKe 


A1=2.-*PO(40)-PD( 38) 
EPS=0.1 

A2=PD(42)-Al 

IF (A2eGTeEPS)IGO TO 70 
IF(DABS(A2)6GTeEPSIGO TO 80 
GG TO 90 

HO=HO-1.0F£-07 

ITER-=ITER+4+1 

TF CITER-100) 555,55, 64 
WRITE(6,100)HO,ITER 

STOP 

HO=HO+1 »~20F-08 

ITER=IFTER+41 
IFCITER-100)55,55,61 
WRITE(6,100)HO,ITER 

STOP 

WRITE(6,100)HO.S ITER 
NO 110 1=1,101 
PD(41)=(PD(42)+4+PD(40))72- 
WRITE(6.120)PD(1) 
FORMAT(T102E£1526) 
CONT ENVE 

ITER=1 
DO 190 J=1,101 

HTC I)=HDO( 1) 
GALS ECALTRs EC eHOsAsPDeHTT151015101499sHD,1) 
CALL CULSHIF{HTT1+.AeR HO) 
CALL ECAL( Rs ECs HGs As PDsHTT 619015101 299sHDs1) 
CALL CLSHIFCHTTs:A,R5HO) 

90 196 1=1,101 

WRITE( 6,197) HTTIC(I)VsHTT(LI),1 
FORMAT(T103s2F1528 »sT60s13) 
CONTINUE 


COMPARE DEFLECTIONS DUE TQ MOMENTUM EQNSeAND ELAST 
* “IC EGNS 


SOLVE "INVERSE "EQUATIONS 
ACONS=PS¥#HDO*®HD*¥2 0006 /(120*A*0 0 4¥1 245031 0F-05 ) 
BCONS=U/2. 
CCONS=U/2.6 
DCONS=12+(02009*PDN(51)/( 126420200 *PD(51) » ) 
DO 130 1=41,61 
PNDC(I)=PO0(1)7PS 
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140 


160 


250 
150 
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1900 
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DO 140 1=41,40 
DPXND CT I=(PNDO(T4+1)—-PND(1-1)9/70e.2 

DPXND(61)=0.64 

J=1 

M=3 

MM=41 

NN=61 

DQ 150 I=MM,NN 

DCT )=1 +10 -009*PD(1)7(1.4+02.0250*PD(1)) ) 
V{(T)=DEXP(PND(1)) 

XCOF (JI)=DCONS*CCGNS 

XCOF (341 )9=-D( 1) *CCONS 

XCOF (J#+2)=0- 

XCOF (53 +3)=ACONS*DPXND(1I)*¥ OCT) /SV(1) 

CALL DPOLRT(XCOF,COF. 3sROOTR,ROOTI,IER) 
WRITE (62160) ROOTR(1).,ROOTI(1).R00TR(2),ROOTI(2).R0 
* OTR (3).R00TI( 3), 
ATER 

FORMAT(6E15-82T100.,13) 

HT CT Y=RGOTRC 1) *HG 

WRITEC6s250)HT(I1) 

FORMAT(T10.E15.8) 

CONTINUE 

HT(61)=HO 

CALL HELMATIHT.ECsP,A) 

CALL ECALICHTT»sHT,PD.ECSA) 


IFLAG=1 

IFLAG=5 

IFC IFLAG.EQe1)GO0 TO 900 

WE WRITE VALUES SG FAR CALCULATED AS DATA IN FILE 
* BRIGU -3 

WRITE (3910)PL5S515S2,E1 2E29R1-2RP 
WRITEC 3,11) ECRs AsPMAX,PS 
FUORMAT(T5.5£61528) 
WRITE(3s100)HO ,ITER 

0G 191 1=1,101 
WRITE(3,192)PD(1),HOCI)SHTT{1I) 
FORMAT (T7105 3E5E1548) 

CONTINUE 

DO 194 17=41,61 

WRITE( 3,193) HT(1) 
FORMAT(T10.,E15.238) 

CONTINUE 

GO TO 901 : 
WE READ BACK DATA FROM FILE 3 
READ(3410)PL59S1,S2sE1»sE22R1sR2 
FORMAT (F10¢592F5 0 352F10 02 s2F5~ 3) 
READ(3011)ECsRsAsPMAXSPS 

READ{3,100)HO,ITER 
FORMAT(E1526.T40.,13) 
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DO 296 T=1,101 

READ 3A¢192)PO0C1) HDC T)IsHTTCI) 

CONTINUE 

DQ 297 1=41,61 

READ( 3s193)HTC(I) 

CONTINUE 

DO 151 I=1,.191 

PDIP{I)=0e 

D(S51)=1200004(0. 009*PD(51))/(1064+002026*PD(51)) 
0G 101 T=1,51 

PSPI(I)=PD(1) 

J=151 

DD 102 1Y=61,.101 

PSPI(II=PD (1) 

J=J+1 

{J=} 

IKN=72 

JJIJ=1 

SPIXK=0.760001D00 

SPI X=0.-800001D00 

Li=25 

LL=35 

il=21 

DO 400 I1G=1,JJjJ 

SPIX=SP1IxX-0.201 

LL=LL+1 

IKN=IKN+4+1 

L=SP1X/0.01 

LPSL=L+51 

LPS2=L4+52 

THE SPICAL ROUTINE FIXES SPIKE POSITION 
0Q 253 1=LP52,151 

PSPIC(I)=0. 

JUK=3 

DO 399 Ni=1.JUK 

CALL SPICAL (Re ECsAsHT sHTTsPD»sD »sPMAX,HO,HSPISPSPI, 
*¥IJ»PDIPsSPIXs»sNI,HSPI1-,IKN) 
CALL DIPCAL CHO, VS»ePSHSPI1 sPSPI»UsAsDeil sPDIP) 
CALL CHEDIP(PDIP,I1J) 

THE DIPCAL ROUTINE FIXES PRESSURE IN DIP 
CONTINUE 

CONTINUE 

RETURN 

END 
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SUBROUTINE SPICAL (Re ECe AsHTeHTT,P0sDsPMAX 2®HO,HSPI 
*sPSPIsIJsPDIP.SPIX NI sHSPIL,IKN) 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION HTT(101)eHT{(101),PD(101)290(101),PSPI(19 
*1)200(0191).HSPI(191) ,HSPI1(191) ePDIP(191)sA1{191), 
*¥PDP(191) 

DIVIDE REGION C AND D INTO 100 EQUAL PARTS 
L=SPIXx 70201 

LPS3=L+53 

LPS1=L+51 

LPS52=L+52 

X=O6 

IFCUNLeGTe1)GO TO 80 

DO 20 1=51,LP51 

Q1l=DABS(1.-X**2) 

PSPICI)=PMAX*¥DSORT(Q1) 

X=X*#029010D000 

READ(S.27)(PSPI(1),1=113.LP53) 

FORMAT(20F 522) 

DO 21 1=1,191 

PDPTIJI=PSPI(1) . 

WE NOW HAVE ASSUMED PRESSURE IN REGION D AS ZERO 


CALCULATE THE ELASTIC DOISPLACEMENT 
CALL ESPCAL(R,ECsAsPSPI2HSPI1,HO) 


DEP CENTRE LINE SHIFT 


DO 75 K=52,151 

DDK 9=1 2 0D00+F( Os CO9*PSPI(K))/(164+0.025*PSPI(K)) 
A1(K)=0(51 ) *HO/DD(K) 

CONTINUE 

C=(CHSP11(052)-HSP11(102))*R/(A*0.50) 

CONS T=HO-HSPI1(52)4C*0.01* ASR 

DO 22 1=52s151 

SHINC=C¥( 1-51 ) *A*0201/R 

HSPI1’1 C1 )=HSPI1C1)4+CGNST4+SHINC 

IN ORDER TO MAKE THE FILM SENSIBLY PARALLEL FILM P 
RESSURES IN REGION C ARE INCREASED 


DO 300 K=52,151 

WRITE (52301) KsPSPI(K) sHSPI1(K),A1(K) 
FORMAT(T19.,135T45.,3F1528) 

CONTINUE 

RETURN 

END 
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SUBROUTINE ESPCAL(RsECsAsPSPI,HSPI,HO) 

IMPLICIT REAL*8{A—H,0-Zz) 

DIMENSION €£(191),U(191),AK(191),PD1(191) »PDD (191 
*)»PSPI(191),HSPI (191 )5SUM2(191)5SUM3(191) 

90 10 J=1,51 

E(J)=(IJ-1)*021-5.e 

DO 15 J=52,4151 

£(J)=(J-51)*0.01 

O8F20.0=152,191 

E(J)=(I-151) 8001 +1. 

=—-5e 

DG 40 K=1,191 

SUM=O0. 

DO 30J=1+189.2 

U( J) =E(9)-x 

U(J*+2)=£(54+2)—xX 

Al=DABS(U{y)) 

A2=1.05F-03 

TF (A1—-A2 125425426 

AK 3} =0% 

eorrTratiz 

AK( J)=025%*U( J) *U( 5) *(DLOG(DABS(U(I)) 9-125) 

M=J+2 

A1=DABS(U(M)) 

TF (A1-A2)18218519 

AK(M)=0- 

GO TO 21 
—AKIM)=005*U(M) *U(M) &( DLOG(DABS(U(M) ) 2-125) 

IF {(J-49)224,22,23 

POL (J )=(—-3 eo *PSPI( I )44.8PSPI( 941 )—-PSPI(9+2))/7( 002) 
PD1(J4#2)=( PSPI(J)-4.*PSPI(J+1)4+3-*PSPI (J4+2) 9/1002 
* ) 

PDD (J)=(PSPI(J)-2-*PSPI(J+1)4PSPI (9429970201 
SUM1L=(PD1 (9) XAK( J)—-PD1 (942) XAK( 942) )-PDD( I) 350 *(U( 
*J) *(AK(J)—-U( I) ®U( 59) 766 D—-U( S42) KC AK (S42)—-UL 42) *UC J 
*+2)/60)) 

GOUTG, 30 

IF (J—149) 24, 24,22 

PD1 (I) =(-36*PSPI(J)44.*PSPI (J41)-PSPI(34+2))7(0202) 
PD1(J+2)=( PSPI(J)—-4.¥*¥PSPI(J41)43.*PSPI(J42) )7(020 
x 2) . 
PDD (I) =(PSPI(J)-2 6 *PSPI( J41)4PSPI(J+2))7020001 
SUML=(PD1{ J) ¥AKOCI)—PD1( IS 42) ¥AK (942) )-—PDD( J )730*(U( 
* J)*CAK(9)—U( J) *U 
AC II/60 Y-U( S42) *( AK ( 542 )—U( S42) U1 S42)/60)) 
SUM=SUM+SUM1 

IF (K-51)41+41,42 

X=X+0e01 

SUM2(K)=SUM 

GO 70°40 

IF (K-151) 43543441 

X=X+0201 

SUM2(K)=SUM 
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CONTINUE 

X=—5e 

DYNO 50 K=1,.51 

SUM3 (K)=SUM2(K)—-SUM2(51) 
GI=(XRKD2 PRC AKKD)D (2. %2) 
G2=2 «6 FECKA¥SUM3(K) 
HSPI(K) =HO+G1—G2 

X=X+0el 

X=Oe6 

DO 60 K=51%151 

SUM3 (K)=SUM2(K )-SUM2(51) 
GIl=CX¥*K2) *(A¥K2)/(2.*R) 
G2=2 e FECXAXSUM3(K) 
HSPI(K)=HO+4+G1—-G2 
X=X+O0e01 

X=1.e 

DU 70 K=151,191 

SUM3 (0K) =SUM2(K)-SUM2(51) 
GI=(X #¥2) €(AK*K 2D (20%R) 
G2=2 e FEC*K¥AKSUM3(K) 
HSPI(K)Y=HO4+G1—-G2 

X=X+0e1 

RETURN 

END 
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SUBROUTINE CHEDIP(POIP,IJ) 
IMPLICIT REAL*8(A-H,0-Z) 
DIMENSTON PDIP(C191).,PTOIPCI91) 
TFC TI-1)10210220 

DO 15 f=1,191 

PTDIPCT)=PDIPCI) 

RETURN 

EPS=1.e0E-02 

SUM=0.6 

DG 30 1=103,150 
AL=DASBS((PTDIPCI)—-PODIPC( IT) /PDIPCT)) 
TF (SUM—Al )36%230.30 

SUM=AI1 

CONTINUE 

ITF (SUM—-EPS)40,40,41 

ITF{ 15-5) 42 543543 

WRITE(6.,44)15 


FORMAT (T5,"*NO CONVERGENCE IN*, 12,5" 1TERATIONS® ) 


GO TG 42 
WRITE(6245)1J5 


FORMAT (TS, *CONVERGED IN *, 12s" ITERATIONS *) 


TJ=5 

00 50 1f=1,.191 
PTDIPCI)=PDIPCI) 
RETURN 

END 
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SUBROUTINE DIPCALCHGs VS2PSeHSPI sPSPILsUsAsDelLLePDIP 
) 

IMPLICIT REAL*8( A—-H,0-Z) 

DIMENSTON HSPI(191).PSPT(191),00101),Z(101)5 VAL(1 


*O1)2HC101).,PDIP(191) 


VAL(1)=0. 
ROO0=D(51) 
UB=6 e ®VS*U*1 2456-057 (PS ¥*A*2000.) 
HOND=HO/A 
RO=1e 
J=153-LL 
H(1)=0(51)*HOSA 
DO 10 f=2etL 
H(ITJ=HSPICLIUISA 
RO=1 2 0D000+(02009*PSPI(J)I/1(1264+0.026*PSPI(9)) 
VAL (CT J=UB* THC T)—HOND*ROO/RDIS CH UI )**3) 
J=J+1 
NDOIM=LL 
E=0.01 
CALL DCSF(E.sVAL»Z,NDIM) 
DO 50 1=1,LtL 
Z(1)=Z(1)41. 
WRITE(6,20)2Z(¢( 1) 
FORMAT(T10.2E15.28) 
CONT INUVE 
J=153-LL 
DO 30 MSH 2.LU 
Al=1le/(1--Z{(1)) 
IF {A1)603,60,70 
A2=DLGG(Al1) 
PDIP(IJ=A2*PS 
WRITE(6s40)A25J5s,PDIP(J) 
FORMAT(T109F15238,T35.,135T55.F1528) 
6D TQ 30 
WRITE(6.,61)LL 
FORMAT(T10%s*NECESSARY TO MOVE SPIKE TG THE RIGHT.V 


* ALUE -OF Li-*.13)} 


STOP 

J=J+1 

J=153-LL 

DO 80 T=J2:150 
PSPI(T)J=POIP(1) 
PSPI(151)=06 

J=152 

PDIP (J-LLI=PSPI(J-LL) 
RETURN 

END 
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SUBROUTINE DPOLRT( XCCF sCOF oMs ROOTR»ROOTI TER) 
IMPLICIT REAL¥*8(A-H, O-Z) 
DIMENSION XCOF (1) sCOF( 1) sROGTR(1) »sROUTI(1) 
IFIT=0 
N=M 
TER=0 
[FCXCOF (N4#1))10225410 
IF(N)15,15,4.32 


SETHRERBGR CODE TD 1 


TER=1 
RETURN 
SET ERROR CODE TO 4 
TER=4 
SORUGEZO 


SETS ERRGRYENDE: TO 2 


TER=2 
GHETS) 20 
IFIN-36) 35%35,30 
NX=M 
NXX=N#1 
N2=1 
KJ1L=N+1 
DQ 40 L=1,KJI1 
MT=KJ1-L+1 
COFUMT)=XCOF{L) 


SET INITIAL VALUES 


XO=042C0500101 
YO=0-01000101 


ZERO INITIAL COUNTER 


IN=0 
X= XO 


INCREMENT INITIAL VALUES AND COUNTER 
XO=-10-20*YO 
YO=-10.20*X 


SET X AND Y TO CURRENT VALUE 


X=XD 
Y=v0 
IN=IN#1 
GO TO 59 
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TS 


105 


IFIT=1 
XPR=X 
YPR=Y 


EVALUATE POLYNOMIAL ANDO DERIVATIVES 
ICT=0 

UxXx=0.0 

UY=0.0 

V=0.0 

YT=0.0 

XT=1.0 

U=CUF (N41) 

IF (U)65.21305365 

DG 70 I=1,N 

L=N-[+1 

TEMP=COF{L) 
XT2=X*XT-VRVT 
YT2=X¥*¥VIT+tY*XT 
U=U4+TEMP*xXT2 
V=V+TEMP#YT2 

FIi=! 

UX=UX4F 1*¥XT* TEMP 
UY=UY-FI*YT*TEMP 
XT=XT2 

Vaan 

SUMS Q=UX ®UX4UY*X UY 

{F (SUMSO)755110.,75 
DX=( VeUY—U*UX) ZSUMSG 
X=X+DX 

DY=—- (U*UY4*V*eUX)I/SUMSO 
Y=Y+DyY 

IF (DABS(DY)+DABS(DX)—-1-00-05)100,80.80 


STEP ITERATICN COUNTER 


ICT=1CT+1 
TFC ICT—500)60 385385 
IF CTIFIT)100590,100 
TFC IN-95)50095395 


SET ERROR CODE TN 3 
IER=3 
GO TO 20 
DCO 105 L=1+NxXx 
MT=KJ1-L+1 
TEMP=XCOF (MT) 
XCOF (MT )=COF(L) 
COF(L)=TEMP 
ITEMP=N 
N=NX 
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110 
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120 
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130 


160 


NX=1 TEMP 
IFC IFIT)1202554120 
IFC LF IT )L152500115 
X=XPR 
Y=YPR 
If IT=0 
TF (DABSCY )-120D-4*DABS{X)) 135.125.125 
ALPHA=X +X 
SUMSQ=X*#X4#Y*Y 
N=N-2 
GO TO 140 
X=0e0 
NX=NX-1 
NXX=NXX—1 
Y=0% 0 
SUMSQ=0.90 
ALPHA=X 
N=N—-1 
COF(2)=COF(2)+ALPHA*COF(1) 
O00 150 L=2.,N 
COF (L#1)=COF(L41)+4AL PHA*COF(L )—SUMSQ*COF(L-1) 
ROOTI(N?2)=Y 
ROOTRIN2 ) =X 
N2=N2+41 
IF (SUMSQ) 160.165.5160 
=-Y 
SUMS G=0.20 
GO TO 155 
TFIN) 20220245 
END 
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SUBROUTINE RRK2CFUNeHse XI oe YI es Je JJ VEC) 


IMPLICIT REAL*8(A—-H.0—-Z) 
DIMENSION VEC(41) 
VEC(1)=0-e 
H2=H/26 
Y=YI 
X=XTI 

DO 2 1=1,jJ 
DO 1 J2=1leJ 
TI=H* FUNCT X.Y) 
T2=HEFUNCXFH2VtTI/26) 
T3=H*FUN(X4H2eVY*T2/2 6) 
T4=H*¥FUN( X4Hs Y+TS) 
Y=Vt( Tit 2c ¥T2+tPe#T3+T4)/60 
X=X+H 
VEC(141)=Y 
RETURN 
END 
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FUNCTIGN FUN{ X,Y) 

IMPLICIT REAL *®8( A-H,0-Z) 

COMMON ReAsHO»vsETA0sPSsPMAX.U 
A1=DOSQRT(OCABS(X**2-1.)) 

HB= (A*A/R) *¥( DABS X¥*A1 )—-DLOG(DABS{I DABS(X)+A1)))/26 
H=HBAHD 

PO=PMAX/SPS 

RBO= 1 et ((2003*PS*POI/S( 1 ote. 0067*#PS#¥PO)) 
HND=H/A 

HBND=HB/A 

HOND=HG/A 

RG=1 e tC (CeOC3EPS*EY)I S( 1Le+20067*PS¥Y)) 

UB=6 6 XE TAN*U* 1245F—-O05/ (PS*A*2000. ) 
FUN=UB*DEXPCY ) ®( HND—HOND* ROO/RO)/ (CHND**3 ) 
RETURN 

END 
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SUBROUTINE ECAL(ResEC sHOsAsPDeHT sKCoKDsKEe HDoL) 
IMPLICIT REAL*8{( A-H,O-Z) 
DIMENSION £01019 .U0101)sAK(101),.PD1(101) »PDD(101), 


*SUM2 (101) .5SUM3(101) sHT(101) sHTT(101) sPD(101),HD(10 
*1) 


WRITE (6,1) REC sHOSA KC eKDe KE aL 
FORMAT(4F15.8.414) 
DQ 10 J=bLskKC 

E(J)=(I-1)*001 —-Se 
CONTINUE 

=-Se 

0G 40 K=L,KD 

SUM=O0-6 

DO 20 J=LeKEs2 
UC JI=ECII-X 
U(I+2)=F(J4+2)-—X 

AL1=DABS(UJ)) 

A2=1.20E-03 

IFCA1—-A2 915915916 

AK{J)=0. 

GO FG 17 
AK(J)=0-5*U( J) *U( J) *(DLOG(DABS(U(J)))-1.-5) 
M=J+2 

A1l=DABS(UIM)) 

IF (A1—-A2)18218519 

AK(M)=0-6 

69 TO 21 
AK(M)=025*U(M) *U(M) *( DLOG(DABS(U(M) ))-1.5) 
PD1(J)=(-3e*PD( I) 442¥*PD( J4l )-PD(IU4F2) 0/1022) 
PDI 542 )=( PD I)-4 0 #*PD( Jt1 +3 6*PD( St2)9/(002) 
POD (JI=(PC(I)-2 -*¥PO( J#1 D4PD(IN42))70201 
SUM1=(PD1 (045) *AK(J)—-PD1 (0542) ¥AK( 942) )-PDD( 9) /3.*(U( 


eI) EC AKC IJ)-U0C I) FUT 60 )-U0 S42) (AK (5 4+2)-U( 942) UL S 
*4+2)/60)) 


SUM=SUM+SUM"1 
[F(KE2EQe299)G0 TO 22 
Al= DABS(E(KC)—X) 
IF (A1—-A2) 22322523 
ADD=0-6 
GO TO 24 
ADD=PO( KC) € (EC KC )—X) ¥(DLOGUCABSTE( KC) —-X) )-1-) 
SUM=SUM+4 ADD 
TF (L-1) 25225526 
SUBTRA =0- 
‘Su B) Nt) 274 
A1L=DABS(E(LI-X) 
TF (A1-A2)28528529 
SUBTRA =06 
GO TO 27 
SUBTRA =PD(L)¥*(E(LI-X) ¥(DLOG(DABS(ETL I-X) )-1.) 
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SUM=SUM-SUBTRA 

KX=X+0e1l 

TF CK eEQe51 )SUMX=SUM 
SUM2 (K)=SUM 

CONT INVE 

X=-5e 

DO SO K=1.KD 
SUM3(K)=SUM2(K)-SUM2(51) 
GL= (XK #2) ¥( AKKS)/( 2e*R) 
G2=2 e KECKAXSUM3BIK) 

HT(K )=HO+G1-G2 

X=X+0e1 

CONTINUF 

RETURN 

END 
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SUBROUTINE CLSHIF CHT sAeReHC) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION HT@(101) 

C= (HT (49)-HT(50) 2¥RI( 0.1 *A) 

CONS T=HO—-HT(50)4(HT (49 )—-HT(50)) 

X=-5e 

DO 195 =1,101 

HTC LD =HHT CI) -CeX¥*¥AYR+ECONST 
195 X=X+0e01 

RETURN 

END 
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SUBROUTINE HELMATCHT»sEC eRe AA) 

IMPLICIT REAL*8( A-K,0-Z) 

DIMENSTON HT (101) eCELPD(1001)2A(10 510) sHCH(1001)58% 
(10)-RCT10.1) 


THE MATRIX FOR STRAIN IN REGION B DUE TO DELTAP IS 
SET 


HH=0e1000 

B(1L)=28e/90¢*DLO0G( 10. *HH) HH 
B(2)=326%*4.*DLOG(9 e*HH) XHH/90. 
8(3)=12.¥*4.*DLOG( 8 e*HH) *HH/90. 
8(4)=32.*4.*DLOG( 7. *HH) #HH/90. 
B(5)=20*7e*4e*DLOG (66 *HH) *HH/90. 

B (6 )=HH*¥ 32 0¥46 *DL0G(50*HH)/90. 

BCT )=HH*120*4.*DLOGi 4.*HH)/90.~ 

8 (8 )=HH*320*4.*DLOG( 3. *HH)/90. 

B(9 )=HH*7 0 *4 6 *DLOG( 2 2 ¥HH)/ 900 tHH* 0 05¥*0L0G( 26 *HH) 
B10 )=HH*0 eS5*DLOG(HH) 
A(1s1)=HH*{DLOG(HH)-1.)-Bi1) 
A(2.1)=-B(1) 
Ai3s1)=0.5*HH*DLOG(2.*HH)-8(1) 
A(4,1)=0+-333333*HH*0LO0G( 36*HH)-8(1) 
A(551)=0.375*HH*D0LOG(46*HH)-B1{1) 

K1=5 

K2=4 

K3=3 

K4=2 

K5=1 

OD) (1,08 =G4 10 
ACIe1)=70%4-e*DLOG( K1*HH) #HH/90--B(1) 
AC1s 2)=320*4e*DLOG(K 2 *HH) *HH/90--B( 2) 
Al1s3)=120*406¥*OLOG(K3*HH) *HH/902-58( 3) 
A( 194 )=322*4 6*DL 0G (K4*HH) *HH/ 906-514) 
Al1+5)=7 0 *4e*DLOG(K5*HH) #HH/90.-B(5) 
K1l=K1+1 

K2=K2+41 

K3=K3+1 

K4=K4+1 
K5=K5+1 

A(7s5)=Al 755) 4005 *HH*OLOG( 26¥*HH) 
A(855)=A( 855) 40-3333 *HH*DLOG (36 *HH) 
A(9¢5)=A(945) +0. 375*HH*DLOG( 4. *HH) 
A(1055)=AC10,5)428e*DLOG(5 e*HHI/90. 

K=2 

DO 15 1=1.5 

AC 1s K)=7 0 #40 *DLOG( HH ) *HH/90 o-B(K) 
Al1sK+1)=320%4 e*HH*¥DLOG( 26#*HH)/90-.—-B(K4t1) 
AC 1s K+2)=120%*4e*DLOG(3 «6 *HH) ¥HH/90 2-B(K+2) 
Al1sK4+3)=320*4e*DLOG( 406 ¥*HH) *HH/90 o—B(K+3) 
ACI 9 K+4)=7e€4e*DLOG( 50 HH) #HH/90.-BI K+4) 
K=K+1 









fa « A \ Pe 5. id avi ful? | 
lA +} 7 7 { Sale 


: j ' if ) TH wO1lee anne: 
(2e01IRe O85) 0 
on 
+ *jyOTae an 
342 * 
1. Ooh 
*? r 11a 


u 


Cidedielnniie on. 


; 41% 66? . SEa1S, FS 3A 
Si steetia 
Pcp sies , 
O94 teeta Pes: 


. i41¢eay 
- 





Seo un=gn 
a? 
item aw ae 
tomNate | OF 
Nhe. SIDQ408HHE OEE ARES MA 
LMn@*.F +o MISMO LEE E were Bbactes ; 
(ivet@. es eit) jeh 4HeSTT.Or le, 


ae 
«DON Cite. 27 20 JG@ abroad ze 
. 7 | = 


i” 


















16 


ig 
20 


H2] 


A(156)=A(156)4282*DLOG(5e*HE) €HH/90} 
A(297)=A(227) +284 *FDLOUG(5 2* HH) ®XHH/90-} 
A(398)=A(398)40037S*HH*EDLOG( Se *HH) 
A(459)=A(44.9 740263333 3 ¥*¥HH*EDLOG(S 0 XHH) 
A(5e010)=A55.10)+00e S*HH*¥OLOG(5 6*HH) 
DG 16 [=2,10 

ACI, 1)=2 6 *HH* DLOGCHH )-HH-BC(T) 


REST OF THE ELEMENTS ARE LISTED INDIVIDUALLY 
A(3s2)=0-e5*HH*OLGG(HH)-B8(2) 
Al4, 2)=1- 3333 3*HH*EDLOG( 2 e*HH)-8B(2) 
A({592)=9e *HH*DLOG( 3e*HH)I/8--B( 2) 
A(4-43)=0233333*HH*DLOG(HH )-B(3) 
A(503)=9 0 *HH*EDLOG( 2 0 #HH)/82-8(3) 
A(544)9=0 e375 *HH*DLCG (HH )-8(4) 
Al7s6)=025*HH*DLOG (HH)-B(6) 
A(8e6 )=1e 33 33%HH¥*¥DLOG{ 2 e ®HHI-B(6) 
A(9¢6)=9 oe *¥HH*DLUG( 30 *HH)/8 2-816) 
A1l1056 )=HH*32 .*4e*DL0G(42*HH)/90.-8(6) 
Alls 79=320*4 e K¥HH*DLOG( He *HHI/902—-B(7) 
A(657)=70*€40e *HH*DLOG(HH)/90.-B(7) 
ABs 7 )=023333*HH*DLOG(HH )-B(7) 
A(99 7 )=9 e X¥HH*FDLOG( 20 #HH) /8--B( 7) 
A(1057)=120%*40*HH¥*¥DLOG(3 e#HH)I/90.—-8(7) 
A(1,8)=120%4e*HH¥EDLOG( 72 *HH)/90.—-8(8) 
Al 2s 8)=326*4.*HH*¥DLOG( 60 ¥HH)/902-8(8) 
Al 6,58 )=320*46*HH*XDLOG( 20 *HH)/902-8( 8) 
A(7+8)=0¢ 375*HH* DLOG (HH )—-B( 8) 
A( 9s 8)=0 0 375*HH*DLOG(HH)-BI(8) 
A(1048)=32 0 *4 ¢e XHH*¥EDLOG( 2 2 FHH)/902-8(3) 
A€1., 9)=322%*4 .*HH*EDLOG{ Be #HHI/90--BI9) 
A(209)=1224 e*®HH*¥DOLOG{ 70 ¥HH)/902-8(9) 
A(3s9)=9e *HH*DLOG(60*HH)/8 --B(9) 
A(629)=12 6*HH*4 o *DL0G( 3. *HH)/90-—-B(9) 
A(7s9)=9e *HH*¥DLOG{ 22*HH) /8 2-819) 
A( 8,9) =HH*DLOG(HH )/3.-8(9) 
A(10s9)=702*4e*HH*OLOG(HH)/902 -8(9) 
A(1510)=HH*7 0 €406¥*¥DLOG( 9 0 HH )/90 0 tHH¥O 0 S¥DLOG( 90 *HH 
»-B(10) 
AC2.10)=32e*40e*HH*DLOG(8 2 *¥HHI/90.—-B110) 
A(3910)=90e *HH*¥DLOG(7 « €HH)/82-8(10) 
A(4510)=40 *HH*¥DLOG(6 oe ®HH)/3--3110) 
Al6.10)=320*4 6 FHH¥EDLOG(4 0 *HH) /902-83(10) 
A17210)=9 6 *HH*¥DLOG( 3 e*#HHI/82-B8(10) 
A(8210)=4 6 *HH*¥DLOG(2 « *HH)/32-B8(10) 
A(9,10)=0.5*HH*DLOG( HH)—-B(10) 
DOVeOY f= 14110 
WRITE (6519) (ACI 95) sJ=1510) 
FORMAT(10F12-8) 
CONT INUE 
K=41 
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100 
25 


D0 25 {=1,10 
RCL»sLIAHTCKIS( 26 *AAKET) 
WRITE(6s1LOOIRGC IT o1). HTE(K) sAASEC 
FORMAT(4F15-8) 

K=K+1 

M=10 

N=1 

EPS=1.0£-03 

CALL DGELG{R»sAMseNsEPSSIER) 
WRITE(6,19)¢% RC1,1)9,1=1,10) 
RETURN 

END 
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SUBROUTINE ECALICHTTeHTePDeEC eA) 

IMPLICIT REAL*3(A-H,P-Z) 

DIMENSION HTT(101) HT CLOLI»sPDOCLOL)®ETLTOLI)sULC LOL) eA 
*K (1101) C1969) .D(9o1)+G (101) 

DO 6 1=41,51 

WRITECEG,5) HTITCI)V SHT(L) 


5 FORMAT(T105,2E15-8) 
a CONTINUE 

WRITE (6e7)YEC 2A 
7 FORMAT(T1092F15.28) 


00 10 J=1.101 

EC J)=( J-1)¥*%021 -Se 
10 CONT TINUE 

X=-121D00 

DG 40 K=159 

X=X+02e1D00 

DO 20 J=41551%2 

U(JI=EL(I)—-X 

U{(J+2)=E(J+2)-xX 

AL=DABS(UJ)) 

A2=1.0£-03 

IF (AL-A2)15,15,16 
15 AK(J)=0.4 

GO TO 17 
16 AK(J)=005*U0I) #U( I) * (OLOG(DABS(U(J)) 9-1-5) 
Wee M=J+2 

A1L=DABS(U(M)) 

IF 1A1-A2)18.13519 

138 AK(M)=0-6 


(Eye) VG) et 
19 AK(M)=0 eS¥*U(M) £U(M)¥*(DOLOG(DABS(U(M)))-12-5) 
21 GUI IV=HUTIV*¥CAK( I) —-UC IS) FUT I) 760 )-U( 942) ¥(AK( 5t+2)-Ul 


* J+2)*U(It2)/60) 
G(JI=GUIIIZ6 
20 CONTINUE 
G(51)=G6(49) 
06 30 1[1=41,51,2 
AK (CI )=AKCT)—-AK(51) 
30 G(1)=G(1I)-G(51 ) 
C1=4 02 FAK (41 9/0 e244 oe FAK 43 ISD e2t2e*®O(41)/0201 
C2=—-AK (41 9/0 262-3 e FAK (43) /022-G(41)/70-01 
C3=—-3 2 *AK(439/0 e2-AK(45)/0 22-G(43)/0.01 
C4=4 6 ¥AK (43 )/0 0244 oe *X¥AK(45)/ 00242 20*G6(43)/0-01 
C5=—-AK (437970 «2-30 *®AK(45)/0 22-G(43)/0-01 
C6=—3e¢*AK(45)/0e2-AK(47)/0 e2-G{(45)/0 e901 
C7=4 @ FAK (45)70 0244 ce ¥AK(47)/0224220*6(45)/0201 
C8=—-AK (45 )/0 02-3 0 X¥AK(47)/0 02-G(45)/0-01 
C9=— 36* AK(47)/022-AK(49) /0022-G6(47)/0-01 
D1=4 6 *AK (47 9/0 0244 ce ¥AK(49 )/ 00242 6%*6(47)/0- 01 
D2=—-AK (47 )/0 22-3 eX AK(49) /0 0 2-G(47)/02-01 
03=—-3 *¥AK(49)/0 2 2-AK(51)/022-G6149) 70-01 


cc 
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40 


$9 
100 


30 


159 


D4=4 6 FAK (49970 0244 eKAK( 51 IS 00242 0 ¥G(49)/0201 


C(K,1)=Cl 
C{Ke2)=C24+C3 
CUK~e3)=C4 
C(Ks4)=C5+CE 
C(UKse SIP =C7 
C(K»6)=C84C2 
C{K,7)=D1 
C(K»8)=D2+4+03 
Ciks9)=D4 
CONT INUE 
J=41 
DG 100 1=1,9 
D(1slLI=CHT TL JI-HT( J) (2 e* ECXA) 
WRITE(6299)D(1I01) 
FORMAT(T10.E1528) 
J=J+1 
OPS=1.0F-03 
M=G 
N=1 
CALL OGELG(D»CesMeNeOPS,IER) 
WRITE(6,90)070( 151) ,T=159) 
FORMAT(T5.9E12-5) 
J=42 
00 150 1=1.9 
POTII=PD(JI-D( 151) 
J=J+1 
RETURN 
END 
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SUBROUTINE THERMA(P,H,TLCI »TLC2) 

THERMAL SOLUTION 

IMPLICIT REAL *8(A-Hs.K—-Z) 

DIMENSION P(200) 6H(200)5A1(200) .81(200),C1(200)5D1 
*(200)%T(200).P1(200).Q1(200) 5G61(200) sDPX(200) sDMU1 
*(200)+YMUBA1( 200) .DMUE( 200) .DMUEP (200), TLC1(200),T 
*LC2( 200) s VL 200) 5 CG(200) » CMEGAN(200) »DOMINT(200) 5 TT 
* (200) sPP(200) »sDDPX(200) sDMUEP2{ 200) »R0B( 200) «UU(290 
*0)sDUY(200).TAU(200) sFTRAC(200 ) sFTRACS(200) «-DFTRAC 
A(200) »DITRAC( 200) 


WRITE (6,9) (PCI) ,T=1,121) 
WRITE(629)CH(1),1T=1,121) 
FORMAT(121F10-.7) 
HO= e 4530-04 
i=1.28370-02 
VIS=0e580-05 
VLS=VIS5 
W=0-32 
EPS2=06 
NONDIMENSTONALISE PRESSURE AND H 


DO 2 1=1,4,121 

PCIJ=P(1T) FL/(22¥*wW) 
H@LI=HC1)/HO 
WRITE(6.9)(PC(1),7T=12121) 
WRITE(5.9)¢(H(1),1=1,121) 
HH=0205 

JDIM=21 

CALL DDET3{(HHePsDPX, JDIM, IER) 
J=1 

DO 3 f{=21,121 

PPCJ)I=P( I) 

J=J+1 

HH=0.005 

JDIM=101 

CALL ODET3(HH»,PP.sDDPX,JDIM,IER) 
DG 4 J=1,101 

OPX( J#20)=DDPX(J) 
WRITE(656,9)(DPX{(1),I=1,121) 


VALUES OF TKE PHYSICAL CONSTANTS ARE READ IN 


READ(5e11)9KSsKL»CPSesCPL,DS,0L 
FORMAT(6F1007) 


CALCULATE VELOCITIES 
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DO 999 I JMK=1.2 
READ (5412) ROLVEL,SLIP 
FORMAT(2F10-5) 
U1I=ROL VEL *(12-SLIPI/S(2-e-SLIP) 
V2=ROLVEL-U1 
WRITECOs1LLOIKS sKL»CPSesCPL»sDSesDLsROLVEL sSLIP»VU2,U1 
FORMATCLOE1L3Z¢4/////) 


CALCULATE REYNOLDS »ePRANDTL»ECKERT NUMBERS 
REY=OS¥ROLVEL¥*¥L/(VIS*32.2%12.- ) 

REYM=REY *HOXxHO/ (L*L) 
PRANL=TPL¥VLS/KL¥4 3200 0% 32 eo 2%*12% 
ECKT=(ROLVEL¥*2)/( CPL ¥S70 oe k€32 0212 0%120¥* 778. ) 
M=WE2O000e S( CPL ¥DL¥*ES700*0 eS ¥*L¥*¥778 ekK126) 
WRITE TE, 1L1LLIIREVYsREYM»PRANL sECKT oM 
FORMAT(5E13-5//7) 

DG 15 T=1,153 

T{(1I)=1-0 

JJ=102 

J1l=1 


TEMPERATURE CALCULATION IN FLUID STREAM BEGINS 


DG 70 JJJ=11,121 
WRITE(6s170)5J5J9 
FORMAT(SOXs* STATION", 13) 

TTER=1 
HRI=KS*HO/KL 
RS1L=L¥DOS*CPS¥EUI FHRI*KHR1¥43200 o/ (LEL*EKS) 
RS2=RS1*¥U2/U1 
WRITE{(6.,112)HRI,RS12RS2 
FORMAT(3E13-5////) 

TF (III-22)20221521 

DXB=0.-05 

GO TO 22 

DXB=0.005 

DYB1=0205 
DYB2=0.205 
DYB=O0 el 

POPS7 “11 it=!.50 
00 25 J1=1.19 

A1(J1 )=3-*RS1442*DXB/(DYBI **2 ) 
RB1l(J1)=-2-*DXB/(DYB1*#2) 
D1i(CJ1)=B1C451) 
C1(I1)V=RS1*(T( SI 451-50) €42-—T(JIt51-101)) 
IF (391 eEQe1 ICICJ1)=C1( J1)-DI(J1) 
CONT INUE 


BOUNDARY CONDITION 
A1L(20)=C1let+ DYB/DYB1) 


Sir 


. : 
| 


a) Se» Jey gomtene) Gasca 

(a9, 0725) TAMHOT 

(ale =. 0 Nit 18 <7 19a ore 
r+ IV eOwe SU 

» 1420 4.199~ 099, See Writ has. 
{eS WONAwE TAOS yTAMAQS 


“ue T25kS, STOWNRI See STARS 
{»S LOS «3h 4 SIV Vv Ja J4v J08820294F 

[se JOM aie YERemyas 

¥¥, 06 SE SOs” ZIV @ 419) 4 MARS 

Se. OTe 999 ta9¢s € jor lJ0R) ep asa 


(.-Si*<ei ® eg (Sy OTAY 1 [277% sOOESORSaN 


M,THIG. MAS, HESS a PIGR Past IA 
(Sin@.) 7300 TARR05 

e412? at) oa 

: et2t241 

CO rect 


<4 tei 


ive wis AE MAT TaAuus j4,3) 2h TARNSRST 


ist, i=404 OF OG 
C4. & (ti «OST tae 
i'r, "HOE TST S*, OST ARNOT 

’ ceraty 


~~. 


rose 2R= LAH 


iA. chee (ane 1 Ghee (UaPegazeese tee. 
PUN Gye Fea ogee 


dng (2H, SO ATE OTT ING. 


CV\AVS E680 P TAMBOR 

its) SeOS.€29-444799 

?0.0=8xa 

*s GT Do 

200 .028K0 | 

*ROLO=s TRTO 

eo, dechva 

 #. 380 

ot, berth 4 

. / Sh _eheit” 
(See (QUOpVanOe. ae (ee ceed 
tiee svar Nex ait 

























St } 





H27 


B1(20)=-1l.e 

D1(20)=- DYB/DYB1 

C1i(20)=0. 

ENERGY EQUATIGN FOR LUBRICANT 

CALL CALMU( ISS »P5TsVLsDMUL »sYMUBAIL »DMUE.,DMUEP.STLCl > 
*¥TLC2 eWel sDMUEP2,U01 ,U2) 

CALL SEEJAYIPRANL »sREYM., DXB, DYB»sHeHOs MsLoEPS2,JI5I/I.W 
*¥ sDPX 5s YMUBA1 »DMUL ~.OMUE eDMVEP 2» GQ eVLS»ROLVEL U1 »U2 VL 
*» OME GAN, OMINT oP SECKT »sROB®UU eDUY»sTAUST) 


on 


On 


30 


40 


41 


4S 


DO 30 J1l=21,29 

A1(J1)=GQ(J1-19) 

BI(J1LI=-DOXB*262 /(DYB*DYBEH( JIS) FHC JIS) ) 
D1i(J1)=B1(0J451) 
CICJLI=ACOMEGANL J1—-19) *®VL{( JI-1L 9) 4D0MINT( J1-19) *T( 59145 


#2) *4e-DMINT(J1-19) *T( YLtl) 


CONT INVE 

BOUNDARY CONDITION BETWEEN LUBRICANT AND SOLID2 
A1l(30)=1.+CYA/SDYB2 

91(30)=-1. 

B1(30)=-DY8&/DY82 

C1(30)=0. 


ENERGY EQUATION FOR SOLID 2 

DO 35 J1=31,49 

A1( 391 )=3e F®RS24+42.¥*¥DOXB/ (DY B82%**2) 
B1(J1)=-2.*DXB/ (DYB2**2) 

Di€@JjJ1)=6810451) 
CACILII=HRS2e(T(II4I1-909%42-T( J J+J1-101)) 
C1(49)=C1(149)-B81(49) 


TRIDIAGONAL MATRIX SOLUTION BEGINS 
P1(1)=A1(1) 
Q1(1)=8101)7P101) 
GI(LI=CI(1)7PL (1) 
DOG 40 J1=2,49 
P10J1)=A1(0J1)-01(051)*Q1(0451-1) 
Bil (J1II=BilCJL)/P1(4I1) 
GICJ1)=(0C1CI1)I-D1 (051) G1 (51-1))7P10I1) 
00 41 i1=1,153 
TTCI)=TC1) 
JJ=102 
J2=49 
J21=J2+1 
TC J214+55)=61 (52) 
DG 45 J1=1,48 
J2=J2-1 
J21=J?+1 
T(J2145I5)=G1( 52 )-O1 (S22) FTC I21+1 455) 
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CHECK FOR CONVERGENCE 


EPS=0.-01 

ERMAX=06 

DOR SSe l= 1%. 3 

QGQ=DABS(C(TCII-TTCIVI/TT(1)) 
Se IF (QQQeGTeERMAX) ERMAX=Q0Q 

IF(CERMAX.GTeEPS)IGO TO 65 

GO TO 69 


CONVERGENCE PARAMETERS ARE PUT IN 


65 DO S56 1=1,153 
56 TUL) =O*eS*TTC(I)0 +6 12-0465) *T( 1) 
ai gd ITER=ITER+1 
69 WRITE(6e171)1TER 
MP An FORMAT(50X, *CONVERGED IN*®',I3s*ITERATIONS® ) 
WRITE(65173)1(T(1),1=123,133) 
Lee3 FORMAT (10Xs*LUBRICANT TEMPERATURES ARE',11F3e5//) 
WRITE (62174)(T(1).1=103.,123) 
174 FORMAT(10X."*LOWER SOLID TEMPERATURES ARE ,21F8e5// 
* ) 
WRITE(60175)(0T(1),1=133,153) 
res. FORMAT(10X,* TOP SOLID TEMPERATURES ARE* .21F825// 
*x =) 
WRITE(6s201)JIJ5J 
201i FORMAT(10X-e*FOR STATION® 613 5%2*%,.T40. *VELUCITY* ,T60 
*&+*DU/DY".T90~2* SHEAR STRESS*,T1055"DENSITY*,T1205°V 
*ISCOSITY*) 
DO 202 IT=1.11 
WRITE(6s203)1,0UU(1I) sDUY( I). TAUCI),ROB(I) »VL( 1) 
203 FORMAT (T5.125T40sF10055T60 oFl 0056 T90sF12 05 eT1055F 1 
*¥ 0242 TI20.E12¢4) 
202 CONT INUE 
WRITE (6,205) 0MUE( JIS) sDMUEP( JIU) sDMUEP2(JJI) s TLC!1 ( 
* JJJ) sTLC2( JIS) 
205 FORMAT(5E12-.5) 
FYTRACCIJSJI=ATAUCLLISH( JIN) 
00 172 {= 1,102 
172 TCT )=T( 1451) 
70 CONT INVE 
DD 176 I=11.1?1 
FYRAC(I-10)=FTRACCI) 
176 CONT INUVE 
HH=0-05 
JOIM=11 
CALL DQSF (HH»sFTRAC sFTRACS, JDIM) 
AX=F TRACS(11) 
DO 300 JKS=21.2121 
300 DFTRAC( JK S—20) =F TRAC(IKS) 
JOIM=101 
HH=0 2005 
CALL DGSF (HH»DFTRAC.DITRACs JDIM) 
AY=DITRAC(101) 
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451 
999 


WRITE C6450) AXKsAY 
FORMAT(2E12 04) 
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FT=CAX+AY) ¥L¥ROLVEL*®VIS/(W*2000 2*HO) 


SVEL=VU2-U1 

WRITE(6,301) SVEL.FT 
FORMAT(10X.*FOR A SLIP VELOCITY OF 
*® ONAL TRACTION IS *5F10¢5//) 

WRITE (65451) (DITRACCI)».I=1.101) 
FORMAT(101€10.24) 

CONTINUE 

STOP 

ENO 


"sF10e05e°FRICTI 
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SUBROUTINE BEEJAY(PRANL,REYM,.DXB, DYBeHsHO»sMslsEPS2 


*¥ »>J JJ 0 WsDPXs YMUBAL Ss DOMUL sO NUE »DMUEP»Q0eVLS»ROLVEL»U1 
*yU2 9 Vis CMEGANs DMINT oP sECKT sROBeUUsDUY e TAUST? 


IMPLICIT REAL*®8(A-H, K-Z) 
DIMENSION H(200) »DPX(200) »s YMUBA1(200) 2DMU1(200).DM 


VE (200) sDOMUEP( 200) 5 QQ( 200) 4 VL( 200) s OMEGAN( 200) sDMINT(200 


id 


30 


dS 


40 


*) sDUY( 200) sUU 200) »P( 200) .RO8( 200) » TAUL 200) sT(200) 


AA1L=2 2 *DXB/((DYB*ee20*(H( III) **2)) 
AA2=WECHC IIS) #¥2 ) KC HOE*2 2 *€20000/( VLS*ROLVEL¥L¥*¥0 05% 


=o) 


AA3=U1/ROLVEL 

UU(1)=AA3 

UU(11)=L2/R0LVEL 

AALO=P(JJIJI) ¥W/(1025¥*L ) 

YY=0.-1D00 

D0 10 1=2,510 

AAS=AA2]*DOPX( ISI) *¥YMUBAICI) 

AAS=(U2-U1)/ROLVEL*( DMUE(IJJI)) ®OMULLI) 

AAE=AA2/ 2 eo *¥DPX( ISI) XOMUE (JSS) /DMUEP( JIS) KOMU1(I) 
UU CT) =AA44+AA5—AABTAA3 

AA7=YY-OMUETIIJI/02¢ ®¥DMUEP(JJJ)) 

DUY ( LJ=AA2Z*DPX( II IIVEAAT/VL (1) 4+(U2-U1 ) SROLVEL*DMUE ( 


* JIJJI/SVL(T) 


YY=YY+021D00 

CONT INUE 

J=123 

DO 30 1=1.511 

ROB (T)=1 +00 009KAALOS( 1040 0026*AA10 )—-0200035*57 00% 


woe (2) — fs J 


J=JI+1 

J=123 

DO 15 I=2,10 

EPS2=0 e00035*57026*(T(J)-1.2) 
J=J+1 

QQCI )=2 « ¥AA1 +36 *¥PRANL*RE YM#ROB( 1) *UU(CT )-2.¥*DXB*HOX 


* HO*¥MXEPS2*UU( I) x 
ADPX( JIJJ) ¥PRANL/{(L*L) 


OMEGANTT )=2 « *DXR*PRANL *ECKT*DUY( 1) *OUY( I) /0H( JJ) * 


* 2) 


DOMINTC I )=PRANL*¥REYM*EROB( I) *UULT) 

CONTINUE ; 

AAT=—-OMUE (IJ II/ (02 © ¥OMUEP( JJ) ) 

DUY(1 J=AA2KDPX( IJJ)*¥AAT/VL(1 9 #(U2-U1) /ROLVEL*OMU 


* ECJIJJISVL(1 ) 


DUYCL J=AAZ¥DPX( JIJS)*AAT/VL(I1 9 40U2-U1)/ROLVEL 
DUY (11 J=AAZKDPX{( ISI) KAAT/VL(11940U2-U1) SROLVEL*DMU 


* EC ISJJII/VL(11) 


DO 40 T=1,11 
TAUCT)=VLC 1) *€0U0Y( 1) 
RETURN 

END 
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SUBROUTINE CALMU(JJJeP oT »eViLeDMU1 s YMUBAL » OMUE » OMUEP 
* »TLCL»sTLC2 oWel sDOMUEP2 oU1 U2) 

IMPLICIT REAL*8(A—-H,K—-Z) 

DIMENSION P(200) 0 T( 200), VL(200) »OMUL (200) »-OMU(200 ) 
* » DMUE (200) »eYMUBA (200) e YMUBAI( 200) »« YMUBA2(200) » YMUB 
*A3(200),TLC1(200)». TLC2(200 ) eDMVEP( 200 ) »DMUEP2 (200) 

AL PHA=—-0 e214*2 o *¥ W/L 

BETA=83006 

GAMA=350e¢*2e* W/L 

AL=ALPHA#P(JJIJ) 

A2=BETA/S70~¢ 

J=123 

090 10 Y=1.11 

AZ=BETA/(T(J)*570.-~) 

A4=GAMA#P (JJIJIIS°I T(J) *57002) 

VLC I =DEXP (AL +A3—A24+A4) 

DMUCT)=1 -/VL( 1) 

J=J+1 

10 CONTINUE 

HH=0.-1D00 

JOIM=11 

CALL DGSF(HHsDMU.DMUL »JDIM) 

DMUE (JJJI=HALe/SOMUICLII) 

YY=06 

900 15 I=1,11 

YMUBACTI=YY/SVLCUI)D 

YMUBA2CII=VYY*¥YY/SVL (1) 

15 YY=YY#0-e1D00 

CALL DOQSF( HH. YMUBA sYMUBA1»JDIM) 

CALL DGSF( HH, YMUBA2.YMUBA3.JOIM) 

OMUEP ( JIS) =H107(20e* YMUBAI(11)) 

DMVEP2(5 JJ I=1 0 /( 30 FYMUBAZ(11)) 

VARAD=DMUE( JJ J) /(3 e XOMUEP2 (J JJ) ) 

SRINU=DMUE(IJJI/Z (02 ¢*¥ OMUEP( JJJ)) 

SRINI=SRINUK*2 

TLC (ISJID HL 20 F(VARAD-SRINI) 

TLO2( JI D=1 0 ¢(U2—-U1 D/C U24U1 ) 01 --D MVE (ISI JI SOMUEPTS 
* JJ)) 

RE TURN 

END 


STOP 0 
08247235 106484 RC=0 
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